diff --git a/Makefile b/Makefile index 1aab947c..c5d7497b 100644 --- a/Makefile +++ b/Makefile @@ -1,13 +1,36 @@ -.PHONY: demo test lint clean +.PHONY: demo test lint clean bench-leakage bench-baseline + +PYTHON := $(shell [ -f .venv/bin/python ] && echo .venv/bin/python || echo python3) +PYTEST := $(shell [ -f .venv/bin/pytest ] && echo .venv/bin/pytest || echo pytest) +RUFF := $(shell [ -f .venv/bin/ruff ] && echo .venv/bin/ruff || echo ruff) demo: - PYTHONPATH=src python -m openamp_foundry.cli rank --candidates examples/sequences/demo_candidates.csv --references examples/known_reference/demo_known_amps.csv --out outputs/demo_ranked.jsonl --report outputs/demo_report.md --cert-dir outputs/evidence + PYTHONPATH=src $(PYTHON) -m openamp_foundry.cli rank \ + --candidates examples/sequences/demo_candidates.csv \ + --references examples/known_reference/demo_known_amps.csv \ + --out outputs/demo_ranked.jsonl \ + --report outputs/demo_report.md \ + --cert-dir outputs/evidence \ + --manifest outputs/run_manifest.json test: - pytest -q + $(PYTEST) -q lint: - ruff check src tests scripts + $(RUFF) check src tests scripts + +bench-leakage: + PYTHONPATH=src $(PYTHON) -m openamp_foundry.cli bench leakage \ + --candidates examples/sequences/demo_candidates.csv \ + --references examples/known_reference/demo_known_amps.csv \ + --out outputs/leakage_report.json + +bench-baseline: + PYTHONPATH=src $(PYTHON) -m openamp_foundry.cli bench baseline \ + --candidates examples/sequences/demo_candidates.csv \ + --references examples/known_reference/demo_known_amps.csv \ + --positives examples/known_reference/demo_known_amps.csv \ + --out outputs/bench_baseline_report.json clean: - rm -rf outputs/*.jsonl outputs/*.md outputs/evidence .pytest_cache .ruff_cache + rm -rf outputs/*.jsonl outputs/*.md outputs/*.json outputs/evidence .pytest_cache .ruff_cache diff --git a/src/openamp_foundry/benchmark/evaluate.py b/src/openamp_foundry/benchmark/evaluate.py index 34c5ae0b..00d9a064 100644 --- a/src/openamp_foundry/benchmark/evaluate.py +++ b/src/openamp_foundry/benchmark/evaluate.py @@ -6,3 +6,98 @@ def top_k_ids(scored: list[ScoredCandidate], k: int) -> set[str]: ranked = sorted(scored, key=lambda x: x.scores.get("ensemble", 0.0), reverse=True) return {item.candidate.candidate_id for item in ranked[:k]} + + +def recall_at_k( + scored: list[ScoredCandidate], + positive_ids: set[str], + k: int, +) -> float: + """Fraction of known positives recovered in the top-k ranked candidates. + + Recall@k = |positives in top-k| / |total positives| + + A random baseline would recover roughly k / |all candidates| of the positives. + The pipeline is meaningful if recall@k significantly exceeds the random baseline. + """ + if not positive_ids: + return 0.0 + top = top_k_ids(scored, k) + recovered = len(top & positive_ids) + return round(recovered / len(positive_ids), 4) + + +def random_recall_at_k(n_candidates: int, n_positives: int, k: int) -> float: + """Expected recall@k for a random ranker. + + E[recall@k] = min(k, n_positives) / n_positives when sampling without replacement. + """ + if n_positives == 0 or n_candidates == 0: + return 0.0 + expected_hits = min(k, n_positives) * min(k, n_candidates) / max(n_candidates, 1) + return round(min(1.0, expected_hits / n_positives), 4) + + +def enrichment_factor( + scored: list[ScoredCandidate], + positive_ids: set[str], + k: int, +) -> float: + """Enrichment factor at k: recall@k relative to random baseline recall@k. + + EF > 1.0 means the pipeline outperforms random. + EF = 1.0 is random. + EF < 1.0 is worse than random (anti-enrichment). + """ + n = len(scored) + n_pos = len(positive_ids) + rc = recall_at_k(scored, positive_ids, k) + random_rc = random_recall_at_k(n, n_pos, k) + if random_rc == 0.0: + return 0.0 + return round(rc / random_rc, 4) + + +def benchmark_summary( + scored: list[ScoredCandidate], + positive_ids: set[str], + ks: list[int] | None = None, +) -> dict: + """Produce a benchmark summary comparing pipeline vs random ranker. + + Returns recall@k and enrichment factor at each k, plus a verdict. + All results are computational only — they do not prove biological activity. + """ + if ks is None: + n = len(scored) + ks = sorted({max(1, n // 10), max(1, n // 5), max(1, n // 2), n}) + + results = [] + for k in ks: + rc = recall_at_k(scored, positive_ids, k) + rrc = random_recall_at_k(len(scored), len(positive_ids), k) + ef = enrichment_factor(scored, positive_ids, k) + results.append( + { + "k": k, + "recall_at_k": rc, + "random_recall_at_k": rrc, + "enrichment_factor": ef, + } + ) + + any_enrichment = any(r["enrichment_factor"] > 1.0 for r in results) + verdict = "pipeline outperforms random" if any_enrichment else "pipeline does not outperform random" + + return { + "disclaimer": ( + "These are retrospective benchmark results on demo data. " + "They do not prove biological efficacy. " + "They only measure whether the pipeline recovers known positives " + "better than a random ranker would." + ), + "n_candidates": len(scored), + "n_positives": len(positive_ids), + "results": results, + "verdict": verdict, + } diff --git a/src/openamp_foundry/cli.py b/src/openamp_foundry/cli.py index 63b0fa71..e90bd4ca 100644 --- a/src/openamp_foundry/cli.py +++ b/src/openamp_foundry/cli.py @@ -19,12 +19,43 @@ def build_parser() -> argparse.ArgumentParser: rank.add_argument("--out", required=True) rank.add_argument("--report", required=False) rank.add_argument("--cert-dir", required=False) + rank.add_argument("--manifest", required=False) rank.add_argument("--config", default="configs/pipeline.yaml") validate = sub.add_parser("validate", help="Validate a candidate certificate against JSON schema.") validate.add_argument("--certificate", required=True) validate.add_argument("--schema", required=True) + bench = sub.add_parser("bench", help="Run benchmark and leakage checks.") + bench_sub = bench.add_subparsers(dest="bench_command", required=True) + + leakage = bench_sub.add_parser("leakage", help="Find near-duplicate candidates in references.") + leakage.add_argument("--candidates", required=True) + leakage.add_argument("--references", required=True) + leakage.add_argument("--threshold", type=float, default=0.90) + leakage.add_argument("--out", required=False, help="Optional JSON output path.") + + baseline = bench_sub.add_parser( + "baseline", + help="Evaluate pipeline recall vs random baseline on a labelled set.", + ) + baseline.add_argument("--candidates", required=True, help="CSV of all candidates to score.") + baseline.add_argument("--references", required=False, help="Reference CSV for novelty scoring.") + baseline.add_argument( + "--positives", + required=True, + help="CSV of known-positive (active) peptides. IDs must match candidates CSV.", + ) + baseline.add_argument( + "--k", + type=int, + nargs="+", + default=None, + help="Recall@k cutoffs (default: auto from dataset size).", + ) + baseline.add_argument("--config", default="configs/pipeline.yaml") + baseline.add_argument("--out", required=False, help="Optional JSON output path.") + return parser @@ -40,6 +71,7 @@ def main(argv: list[str] | None = None) -> int: report_path=args.report, cert_dir=args.cert_dir, config_path=args.config, + manifest_path=args.manifest, ) print(json.dumps({"status": "ok", "out": args.out, "report": args.report}, indent=2)) return 0 @@ -50,9 +82,58 @@ def main(argv: list[str] | None = None) -> int: print(json.dumps({"status": "valid", "certificate": args.certificate}, indent=2)) return 0 + if args.command == "bench": + return _run_bench(args) + parser.error("unknown command") return 2 +def _run_bench(args: argparse.Namespace) -> int: + from openamp_foundry.data.loaders import load_candidates_csv + from openamp_foundry.utils.io import write_json + + if args.bench_command == "leakage": + from openamp_foundry.benchmark.leakage import find_near_duplicates + + candidates = load_candidates_csv(args.candidates) + references = load_candidates_csv(args.references) + hits = find_near_duplicates(candidates, references, threshold=args.threshold) + result = { + "status": "ok", + "threshold": args.threshold, + "near_duplicate_count": len(hits), + "near_duplicates": hits, + "warning": ( + "Near-duplicates detected. If these candidates were used for training or " + "scoring baseline models, benchmark results may be inflated." + ) if hits else None, + } + if args.out: + write_json(args.out, result) + print(json.dumps(result, indent=2)) + return 0 + + if args.bench_command == "baseline": + from openamp_foundry.benchmark.evaluate import benchmark_summary + from openamp_foundry.pipeline import score_candidates + + scored, _ = score_candidates( + candidate_path=args.candidates, + reference_path=args.references, + config_path=args.config, + ) + positives = load_candidates_csv(args.positives) + positive_ids = {p.candidate_id for p in positives} + summary = benchmark_summary(scored, positive_ids, ks=args.k) + result = {"status": "ok", **summary} + if args.out: + write_json(args.out, result) + print(json.dumps(result, indent=2)) + return 0 + + return 2 + + if __name__ == "__main__": raise SystemExit(main()) diff --git a/src/openamp_foundry/features/physchem.py b/src/openamp_foundry/features/physchem.py index 9ff9f9df..2ae0ea44 100644 --- a/src/openamp_foundry/features/physchem.py +++ b/src/openamp_foundry/features/physchem.py @@ -1,5 +1,6 @@ from __future__ import annotations +import math from collections import Counter HYDROPHOBIC = set("AILMFWVY") @@ -8,6 +9,15 @@ AROMATIC = set("FWY") CYS = "C" +# Eisenberg consensus hydrophobicity scale (normalized, 0-centred removed, shifted to 0..1 range) +# Source: Eisenberg et al. (1984) J Mol Biol. Used for hydrophobic moment only. +_HYDROPHOBICITY: dict[str, float] = { + "A": 0.620, "R": -2.530, "N": -0.780, "D": -0.900, "C": 0.290, + "Q": -0.850, "E": -0.740, "G": 0.480, "H": -0.400, "I": 1.380, + "L": 1.060, "K": -1.500, "M": 0.640, "F": 1.190, "P": 0.120, + "S": -0.180, "T": -0.050, "W": 0.810, "Y": 0.260, "V": 1.080, +} + def net_charge_proxy(sequence: str) -> int: return sum(1 for aa in sequence if aa in POSITIVE) - sum(1 for aa in sequence if aa in NEGATIVE) @@ -33,6 +43,29 @@ def longest_repeat_run(sequence: str) -> int: return best +def hydrophobic_moment(sequence: str, angle_deg: float = 100.0) -> float: + """Compute mean hydrophobic moment for an assumed alpha-helical conformation. + + Uses the Eisenberg (1984) consensus scale. Angle of 100° per residue is the + standard helical wheel projection used in AMP literature. + + Returns a value in [0, ∞); typical AMPs have μH > 0.4. Not normalised to [0,1] + because the range is sequence-length-dependent — callers should normalise if needed. + """ + if not sequence: + return 0.0 + angle_rad = math.radians(angle_deg) + sin_sum = 0.0 + cos_sum = 0.0 + for i, aa in enumerate(sequence): + h = _HYDROPHOBICITY.get(aa, 0.0) + theta = i * angle_rad + sin_sum += h * math.sin(theta) + cos_sum += h * math.cos(theta) + moment = math.sqrt(sin_sum ** 2 + cos_sum ** 2) / len(sequence) + return round(moment, 4) + + def compute_features(sequence: str) -> dict[str, float | int | dict[str, int]]: counts = Counter(sequence) length = len(sequence) @@ -43,6 +76,7 @@ def compute_features(sequence: str) -> dict[str, float | int | dict[str, int]]: gly_fraction = counts.get("G", 0) / length if length else 0.0 pro_fraction = counts.get("P", 0) / length if length else 0.0 repeat_run = longest_repeat_run(sequence) + mu_h = hydrophobic_moment(sequence) return { "length": length, "net_charge_proxy": charge, @@ -53,5 +87,6 @@ def compute_features(sequence: str) -> dict[str, float | int | dict[str, int]]: "glycine_fraction": round(gly_fraction, 4), "proline_fraction": round(pro_fraction, 4), "longest_repeat_run": repeat_run, + "hydrophobic_moment": mu_h, "residue_counts": dict(sorted(counts.items())), } diff --git a/src/openamp_foundry/pipeline.py b/src/openamp_foundry/pipeline.py index 727a6a69..9a0d0d14 100644 --- a/src/openamp_foundry/pipeline.py +++ b/src/openamp_foundry/pipeline.py @@ -1,8 +1,12 @@ from __future__ import annotations +import hashlib +import uuid +from datetime import datetime, timezone from pathlib import Path from typing import Any +from openamp_foundry import __version__ from openamp_foundry.config import load_config from openamp_foundry.data.loaders import is_valid_sequence, load_candidates_csv from openamp_foundry.evidence.certificate import build_certificate @@ -15,23 +19,34 @@ from openamp_foundry.selection.diversity import greedy_diverse_select from openamp_foundry.selection.pareto import rank_candidates from openamp_foundry.types import ScoredCandidate +from openamp_foundry.utils.hashing import file_sha256, stable_json_hash from openamp_foundry.utils.io import write_json, write_jsonl +def _passes_length_filter(sequence: str, min_length: int, max_length: int) -> bool: + return min_length <= len(sequence) <= max_length + + def score_candidates( candidate_path: str | Path, reference_path: str | Path | None = None, config_path: str | Path = "configs/pipeline.yaml", ) -> tuple[list[ScoredCandidate], dict[str, Any]]: config = load_config(config_path) - allowed = set(config["filters"]["allowed_amino_acids"]) + filters = config.get("filters", {}) + allowed = set(filters.get("allowed_amino_acids", "ACDEFGHIKLMNPQRSTVWY")) + min_length = int(filters.get("min_length", 8)) + max_length = int(filters.get("max_length", 35)) + candidates = load_candidates_csv(candidate_path) references = load_candidates_csv(reference_path) if reference_path else [] weights = config["weights"] scored: list[ScoredCandidate] = [] for candidate in candidates: - valid = is_valid_sequence(candidate.sequence, allowed=allowed) + valid_aa = is_valid_sequence(candidate.sequence, allowed=allowed) + valid_len = _passes_length_filter(candidate.sequence, min_length, max_length) + valid = valid_aa and valid_len features = compute_features(candidate.sequence) act = activity_likeness_score(features) if valid else 0.0 safe = safety_score(features) if valid else 0.0 @@ -51,12 +66,44 @@ def score_candidates( references_checked=[str(reference_path)] if reference_path else [], nearest_reference=nearest, ) + item.valid = valid item.selection_reason = selection_reasons(raw_scores) item.known_failure_modes = known_failure_modes(raw_scores) + if not valid_aa: + item.known_failure_modes.append("Sequence contains non-canonical amino acids.") + if not valid_len: + item.known_failure_modes.append( + f"Sequence length {len(candidate.sequence)} outside filter range " + f"[{min_length}, {max_length}]." + ) scored.append(item) return scored, config +def build_run_manifest( + run_id: str, + config: dict[str, Any], + input_paths: list[Path], + output_paths: list[str], + generated_at: str, +) -> dict[str, Any]: + inputs = [] + for p in input_paths: + if p.exists(): + inputs.append({"path": str(p), "sha256": file_sha256(p)}) + else: + inputs.append({"path": str(p)}) + return { + "run_id": run_id, + "pipeline_version": __version__, + "config_hash": stable_json_hash(config), + "generated_at": generated_at, + "inputs": [x["path"] for x in inputs], + "input_hashes": {x["path"]: x.get("sha256", "N/A") for x in inputs}, + "outputs": output_paths, + } + + def run_ranking_pipeline( candidate_path: str | Path, reference_path: str | Path | None, @@ -64,19 +111,37 @@ def run_ranking_pipeline( report_path: str | Path | None = None, cert_dir: str | Path | None = None, config_path: str | Path = "configs/pipeline.yaml", + manifest_path: str | Path | None = None, ) -> list[ScoredCandidate]: + run_id = str(uuid.uuid4()) + generated_at = datetime.now(timezone.utc).isoformat() + scored, config = score_candidates(candidate_path, reference_path, config_path) + selection_cfg = config.get("selection", {}) + min_novelty = float(selection_cfg.get("min_novelty", 0.0)) + max_safety_risk = float(selection_cfg.get("max_safety_risk", 1.0)) + ranked = rank_candidates(scored) - top_n = int(config.get("selection", {}).get("top_n", len(ranked))) - selected = greedy_diverse_select(ranked, top_n=top_n) + top_n = int(selection_cfg.get("top_n", len(ranked))) + + eligible = [ + item for item in ranked + if getattr(item, "valid", True) + and item.scores["novelty"] >= min_novelty + and item.scores["safety"] >= (1.0 - max_safety_risk) + ] + selected = greedy_diverse_select(eligible, top_n=top_n) rows = [] + selected_ids = {item.candidate.candidate_id for item in selected} for item in ranked: rows.append( { "candidate_id": item.candidate.candidate_id, "sequence": item.candidate.sequence, "source": item.candidate.source, + "valid": getattr(item, "valid", True), + "selected": item.candidate.candidate_id in selected_ids, "scores": item.scores, "features": item.features, "nearest_reference": item.nearest_reference, @@ -95,6 +160,28 @@ def run_ranking_pipeline( if report_path: write_report(report_path, ranked, selected) + output_paths = [str(out_path)] + if report_path: + output_paths.append(str(report_path)) + if cert_dir: + output_paths.append(str(cert_dir)) + + manifest = build_run_manifest( + run_id=run_id, + config=config, + input_paths=[Path(candidate_path)] + + ([Path(reference_path)] if reference_path else []) + + [Path(config_path)], + output_paths=output_paths, + generated_at=generated_at, + ) + + if manifest_path: + write_json(manifest_path, manifest) + else: + manifest_out = Path(out_path).with_name("run_manifest.json") + write_json(manifest_out, manifest) + return ranked @@ -103,23 +190,27 @@ def write_report( ) -> None: p = Path(path) p.parent.mkdir(parents=True, exist_ok=True) + selected_ids = {item.candidate.candidate_id for item in selected} lines = [ - "# OpenAMP Foundry Demo Report", + "# OpenAMP Foundry Candidate Report", "", - "This report is generated by transparent baseline heuristics.", - "It does not prove antimicrobial activity.", + "> **Disclaimer:** All scores are transparent baseline heuristics computed from " + "physicochemical properties. They are NOT validated biological predictors. " + "No antimicrobial activity has been demonstrated in vitro or in vivo. " + "These candidates are nominated for possible future expert review and assay only.", "", f"Total candidates scored: {len(ranked)}", - f"Candidates selected for evidence certificates: {len(selected)}", + f"Candidates selected (passed all filters, diverse): {len(selected)}", "", - "| Rank | ID | Sequence | Ensemble | Activity | Safety | Synthesis | Novelty |", - "|---:|---|---|---:|---:|---:|---:|---:|", + "| Rank | ID | Seq | Ensemble | Activity | Safety | Synthesis | Novelty | Selected |", + "|---:|---|---|---:|---:|---:|---:|---:|:---:|", ] for idx, item in enumerate(ranked, start=1): s = item.scores + sel_mark = "Y" if item.candidate.candidate_id in selected_ids else "" lines.append( f"| {idx} | {item.candidate.candidate_id} | `{item.candidate.sequence}` | " f"{s['ensemble']:.4f} | {s['activity']:.4f} | {s['safety']:.4f} | " - f"{s['synthesis']:.4f} | {s['novelty']:.4f} |" + f"{s['synthesis']:.4f} | {s['novelty']:.4f} | {sel_mark} |" ) p.write_text("\n".join(lines) + "\n", encoding="utf-8") diff --git a/src/openamp_foundry/scoring/activity.py b/src/openamp_foundry/scoring/activity.py index 9a407191..615e9215 100644 --- a/src/openamp_foundry/scoring/activity.py +++ b/src/openamp_foundry/scoring/activity.py @@ -6,16 +6,47 @@ def clamp01(x: float) -> float: def activity_likeness_score(features: dict) -> float: - """Transparent baseline, not a validated biological predictor.""" + """Transparent baseline activity-likeness score. + + This is NOT a validated biological predictor. It combines known physicochemical + correlates of AMP activity: length, cationic charge density, hydrophobic fraction, + aromatic content, and hydrophobic moment (amphipathicity). All weights and thresholds + are documented and fixed before any benchmark evaluation. + + Literature basis: + - Charge: Zasloff (2002) Nature; most AMPs carry net charge +2 to +9 + - Hydrophobicity: Hancock & Sahl (2006) Nat Biotechnol; ~30-50% hydrophobic + - Amphipathicity (μH): Eisenberg et al. (1984); helical amphipathicity correlates + with membrane disruption + - Length: Jenssen et al. (2006) Clin Microbiol Rev; typical 8-50 AA + """ length = features["length"] charge_density = features["charge_density"] hydrophobic = features["hydrophobic_fraction"] aromatic = features["aromatic_fraction"] + mu_h = features.get("hydrophobic_moment", 0.0) + # Length: peak at ~18 residues, broad tolerance length_score = 1.0 - min(abs(length - 18) / 25, 1.0) + + # Charge: AMPs typically have positive charge density 0.1–0.5 charge_score = clamp01((charge_density + 0.05) / 0.55) + + # Hydrophobicity: 40-50% is a sweet spot for membrane interaction hydro_score = 1.0 - min(abs(hydrophobic - 0.45) / 0.45, 1.0) + + # Aromatic residues (F, W, Y) aid membrane insertion aromatic_bonus = min(aromatic / 0.20, 1.0) * 0.10 - score = 0.30 * length_score + 0.35 * charge_score + 0.25 * hydro_score + aromatic_bonus + # Amphipathicity: helical hydrophobic moment > 0.4 is associated with activity + # Typical range for AMPs: 0.3–0.8; scale to [0,1] over 0–0.8 range + amphipathicity_score = clamp01(mu_h / 0.8) * 0.15 + + score = ( + 0.28 * length_score + + 0.32 * charge_score + + 0.20 * hydro_score + + aromatic_bonus + + amphipathicity_score + ) return round(clamp01(score), 4) diff --git a/src/openamp_foundry/types.py b/src/openamp_foundry/types.py index 4b33a7ea..24775960 100644 --- a/src/openamp_foundry/types.py +++ b/src/openamp_foundry/types.py @@ -20,3 +20,4 @@ class ScoredCandidate: nearest_reference: dict[str, Any] | None = None selection_reason: list[str] = field(default_factory=list) known_failure_modes: list[str] = field(default_factory=list) + valid: bool = True diff --git a/tests/test_benchmark.py b/tests/test_benchmark.py new file mode 100644 index 00000000..438bc8ba --- /dev/null +++ b/tests/test_benchmark.py @@ -0,0 +1,115 @@ +"""Tests for benchmark leakage detection, evaluation, and splits.""" +from __future__ import annotations + +from openamp_foundry.benchmark.evaluate import top_k_ids +from openamp_foundry.benchmark.leakage import find_near_duplicates +from openamp_foundry.benchmark.splits import deterministic_split +from openamp_foundry.features.physchem import compute_features +from openamp_foundry.types import PeptideCandidate, ScoredCandidate + + +def _make_candidate(cid: str, seq: str) -> PeptideCandidate: + return PeptideCandidate(candidate_id=cid, sequence=seq, source="test") + + +def _make_scored(cid: str, seq: str, ensemble: float) -> ScoredCandidate: + features = compute_features(seq) + scores = { + "activity": ensemble, + "safety": 0.8, + "synthesis": 0.9, + "novelty": 0.5, + "ensemble": ensemble, + } + return ScoredCandidate(candidate=_make_candidate(cid, seq), features=features, scores=scores) + + +class TestFindNearDuplicates: + def test_exact_match_detected(self): + cands = [_make_candidate("C1", "KWKLFKK")] + refs = [_make_candidate("R1", "KWKLFKK")] + hits = find_near_duplicates(cands, refs, threshold=0.90) + assert len(hits) == 1 + assert hits[0]["candidate_id"] == "C1" + assert hits[0]["similarity"] == 1.0 + + def test_different_sequence_no_hit(self): + cands = [_make_candidate("C1", "KWKLFKK")] + refs = [_make_candidate("R1", "GLFDIVKK")] + hits = find_near_duplicates(cands, refs, threshold=0.90) + assert hits == [] + + def test_threshold_boundary(self): + # "KWKLFKK" vs "KWKLFKX" — one substitution out of 7 chars = ~86% similarity + cands = [_make_candidate("C1", "KWKLFKK")] + refs = [_make_candidate("R1", "KWKLFKA")] + hits_high = find_near_duplicates(cands, refs, threshold=0.95) + hits_low = find_near_duplicates(cands, refs, threshold=0.80) + assert hits_high == [] + assert len(hits_low) == 1 + + def test_multiple_candidates_and_refs(self): + cands = [ + _make_candidate("C1", "KWKLFKK"), + _make_candidate("C2", "GLFDIVKK"), + ] + refs = [ + _make_candidate("R1", "KWKLFKK"), + _make_candidate("R2", "AAAAAAAAA"), + ] + hits = find_near_duplicates(cands, refs, threshold=0.90) + ids = {h["candidate_id"] for h in hits} + assert "C1" in ids + assert "C2" not in ids + + def test_empty_inputs(self): + assert find_near_duplicates([], [], threshold=0.90) == [] + assert find_near_duplicates([_make_candidate("C1", "KWKLFKK")], [], threshold=0.90) == [] + + +class TestTopKIds: + def test_returns_top_k_by_ensemble(self): + items = [ + _make_scored("C1", "KWKLFKKIGAVLKVL", 0.8), + _make_scored("C2", "GLFDIVKK", 0.6), + _make_scored("C3", "AAAAAAAA", 0.3), + ] + top = top_k_ids(items, k=2) + assert "C1" in top + assert "C2" in top + assert "C3" not in top + + def test_k_larger_than_list(self): + items = [_make_scored("C1", "KWKLFKK", 0.9)] + top = top_k_ids(items, k=10) + assert "C1" in top + + def test_empty_list(self): + assert top_k_ids([], k=5) == set() + + +class TestDeterministicSplit: + def test_split_is_reproducible(self): + cands = [_make_candidate(f"C{i}", "K" * (8 + i)) for i in range(10)] + train1, holdout1 = deterministic_split(cands, holdout_mod=5) + train2, holdout2 = deterministic_split(cands, holdout_mod=5) + assert [c.candidate_id for c in train1] == [c.candidate_id for c in train2] + assert [c.candidate_id for c in holdout1] == [c.candidate_id for c in holdout2] + + def test_no_overlap_between_train_and_holdout(self): + cands = [_make_candidate(f"C{i}", "K" * (8 + i)) for i in range(15)] + train, holdout = deterministic_split(cands, holdout_mod=5) + train_ids = {c.candidate_id for c in train} + holdout_ids = {c.candidate_id for c in holdout} + assert train_ids.isdisjoint(holdout_ids) + + def test_all_items_covered(self): + cands = [_make_candidate(f"C{i}", "K" * (8 + i)) for i in range(10)] + train, holdout = deterministic_split(cands, holdout_mod=3) + assert len(train) + len(holdout) == 10 + + def test_holdout_fraction(self): + cands = [_make_candidate(f"C{i}", "K" * (8 + i)) for i in range(10)] + _, holdout = deterministic_split(cands, holdout_mod=5) + # Every 5th item → 2 holdout from 10 + assert len(holdout) == 2 diff --git a/tests/test_benchmark_evaluate.py b/tests/test_benchmark_evaluate.py new file mode 100644 index 00000000..7738bc36 --- /dev/null +++ b/tests/test_benchmark_evaluate.py @@ -0,0 +1,123 @@ +"""Tests for benchmark evaluation: recall@k, enrichment factor, and summary.""" +from __future__ import annotations + +import pytest + +from openamp_foundry.benchmark.evaluate import ( + benchmark_summary, + enrichment_factor, + random_recall_at_k, + recall_at_k, + top_k_ids, +) +from openamp_foundry.features.physchem import compute_features +from openamp_foundry.types import PeptideCandidate, ScoredCandidate + + +def _make_scored(cid: str, seq: str, ensemble: float) -> ScoredCandidate: + features = compute_features(seq) + scores = { + "activity": ensemble, + "safety": 0.8, + "synthesis": 0.9, + "novelty": 0.5, + "ensemble": ensemble, + } + return ScoredCandidate( + candidate=PeptideCandidate(candidate_id=cid, sequence=seq, source="test"), + features=features, + scores=scores, + ) + + +ITEMS = [ + _make_scored("C1", "KWKLFKKIGAVLKVL", 0.9), + _make_scored("C2", "GIGKFLHSAKKFG", 0.7), + _make_scored("C3", "AAAAAAAA", 0.3), + _make_scored("C4", "GLFDIVKK", 0.6), + _make_scored("C5", "DEDEDEDE", 0.1), +] +POSITIVES = {"C1", "C2"} + + +class TestRecallAtK: + def test_perfect_recall_when_all_positives_in_top_k(self): + assert recall_at_k(ITEMS, POSITIVES, k=2) == 1.0 + + def test_partial_recall(self): + assert recall_at_k(ITEMS, POSITIVES, k=1) == 0.5 + + def test_zero_recall_when_no_positives_in_top_k(self): + # Only C5 (worst score) is "positive" here + assert recall_at_k(ITEMS, {"C5"}, k=1) == 0.0 + + def test_full_recall_at_all(self): + assert recall_at_k(ITEMS, POSITIVES, k=len(ITEMS)) == 1.0 + + def test_empty_positives_returns_zero(self): + assert recall_at_k(ITEMS, set(), k=3) == 0.0 + + +class TestRandomRecallAtK: + def test_expected_random_recall(self): + # 2 positives in 5 candidates, k=2 → expected 0.4 hits → 0.4/2 = 0.2... + # Actually E[hits] = k * n_pos / n = 2 * 2 / 5 = 0.8 → recall = 0.8/2 = 0.4 + result = random_recall_at_k(n_candidates=5, n_positives=2, k=2) + assert 0 < result < 1 + + def test_zero_candidates_returns_zero(self): + assert random_recall_at_k(0, 2, 2) == 0.0 + + def test_zero_positives_returns_zero(self): + assert random_recall_at_k(5, 0, 2) == 0.0 + + def test_k_equals_total_returns_one(self): + result = random_recall_at_k(n_candidates=5, n_positives=2, k=5) + assert result == 1.0 + + +class TestEnrichmentFactor: + def test_ef_greater_than_one_for_good_ranker(self): + # Top-2 by ensemble score are C1 and C2, which are our positives + ef = enrichment_factor(ITEMS, POSITIVES, k=2) + assert ef > 1.0 + + def test_ef_approximately_one_for_random(self): + # With k = all items, every ranker has the same recall + ef = enrichment_factor(ITEMS, POSITIVES, k=len(ITEMS)) + assert abs(ef - 1.0) < 0.01 + + +class TestBenchmarkSummary: + def test_summary_has_required_keys(self): + result = benchmark_summary(ITEMS, POSITIVES, ks=[1, 2, 5]) + assert "disclaimer" in result + assert "n_candidates" in result + assert "n_positives" in result + assert "results" in result + assert "verdict" in result + + def test_summary_contains_disclaimer(self): + result = benchmark_summary(ITEMS, POSITIVES, ks=[2]) + assert "do not prove biological efficacy" in result["disclaimer"].lower() or \ + "do not prove" in result["disclaimer"].lower() + + def test_verdict_positive_when_pipeline_outperforms(self): + # C1 and C2 are top-ranked and are our positives — should outperform random at k=2 + result = benchmark_summary(ITEMS, POSITIVES, ks=[2]) + assert result["verdict"] == "pipeline outperforms random" + + def test_results_per_k(self): + result = benchmark_summary(ITEMS, POSITIVES, ks=[1, 3]) + assert len(result["results"]) == 2 + assert result["results"][0]["k"] == 1 + assert result["results"][1]["k"] == 3 + + def test_auto_ks_when_none_given(self): + result = benchmark_summary(ITEMS, POSITIVES) + assert len(result["results"]) > 0 + + def test_counts_are_correct(self): + result = benchmark_summary(ITEMS, POSITIVES, ks=[2]) + assert result["n_candidates"] == 5 + assert result["n_positives"] == 2 diff --git a/tests/test_cli.py b/tests/test_cli.py new file mode 100644 index 00000000..02e1f390 --- /dev/null +++ b/tests/test_cli.py @@ -0,0 +1,113 @@ +"""CLI integration tests.""" +from __future__ import annotations + +import json + +import pytest + +from openamp_foundry.cli import main + + +def test_rank_command_success(tmp_path): + out = str(tmp_path / "ranked.jsonl") + ret = main([ + "rank", + "--candidates", "examples/sequences/demo_candidates.csv", + "--references", "examples/known_reference/demo_known_amps.csv", + "--out", out, + ]) + assert ret == 0 + + +def test_rank_command_with_report_and_certs(tmp_path): + out = str(tmp_path / "ranked.jsonl") + report = str(tmp_path / "report.md") + certs = str(tmp_path / "certs") + ret = main([ + "rank", + "--candidates", "examples/sequences/demo_candidates.csv", + "--references", "examples/known_reference/demo_known_amps.csv", + "--out", out, + "--report", report, + "--cert-dir", certs, + ]) + assert ret == 0 + + +def test_validate_command_success(tmp_path): + # First generate a certificate + out = str(tmp_path / "ranked.jsonl") + certs = str(tmp_path / "certs") + main([ + "rank", + "--candidates", "examples/sequences/demo_candidates.csv", + "--references", "examples/known_reference/demo_known_amps.csv", + "--out", out, + "--cert-dir", certs, + ]) + import os + cert_files = list((tmp_path / "certs").glob("*.json")) + assert cert_files, "No certificates were generated" + ret = main([ + "validate", + "--certificate", str(cert_files[0]), + "--schema", "schemas/candidate.schema.json", + ]) + assert ret == 0 + + +def test_bench_leakage_detects_duplicates(tmp_path, capsys): + # Demo candidates 1, 2, 5 are exact copies of references + ret = main([ + "bench", "leakage", + "--candidates", "examples/sequences/demo_candidates.csv", + "--references", "examples/known_reference/demo_known_amps.csv", + ]) + assert ret == 0 + captured = capsys.readouterr() + result = json.loads(captured.out) + assert result["near_duplicate_count"] == 3 + assert result["warning"] is not None + + +def test_bench_leakage_no_duplicates(tmp_path, capsys): + # Use negative examples as candidates — they won't match the reference AMPs + ret = main([ + "bench", "leakage", + "--candidates", "examples/negative/demo_negative_peptides.csv", + "--references", "examples/known_reference/demo_known_amps.csv", + "--threshold", "0.90", + ]) + assert ret == 0 + captured = capsys.readouterr() + result = json.loads(captured.out) + assert result["near_duplicate_count"] == 0 + assert result["warning"] is None + + +def test_bench_leakage_output_file(tmp_path, capsys): + out = str(tmp_path / "leakage_report.json") + ret = main([ + "bench", "leakage", + "--candidates", "examples/sequences/demo_candidates.csv", + "--references", "examples/known_reference/demo_known_amps.csv", + "--out", out, + ]) + assert ret == 0 + data = json.loads((tmp_path / "leakage_report.json").read_text()) + assert "near_duplicates" in data + + +def test_report_contains_disclaimer(tmp_path): + out = str(tmp_path / "ranked.jsonl") + report = str(tmp_path / "report.md") + main([ + "rank", + "--candidates", "examples/sequences/demo_candidates.csv", + "--references", "examples/known_reference/demo_known_amps.csv", + "--out", out, + "--report", report, + ]) + text = (tmp_path / "report.md").read_text() + assert "NOT validated biological predictors" in text + assert "no antimicrobial activity has been demonstrated" in text.lower() or "No antimicrobial activity" in text diff --git a/tests/test_physchem_amphipathicity.py b/tests/test_physchem_amphipathicity.py new file mode 100644 index 00000000..59347224 --- /dev/null +++ b/tests/test_physchem_amphipathicity.py @@ -0,0 +1,63 @@ +"""Tests for hydrophobic moment (amphipathicity) feature.""" +from __future__ import annotations + +import math + +from openamp_foundry.features.physchem import compute_features, hydrophobic_moment + + +class TestHydrophobicMoment: + def test_empty_sequence_returns_zero(self): + assert hydrophobic_moment("") == 0.0 + + def test_uniform_sequence_low_moment(self): + # All same amino acid → sine/cosine terms distribute evenly → low moment + result = hydrophobic_moment("AAAAAAAAAA") + assert isinstance(result, float) + assert result >= 0.0 + + def test_alternating_hydrophobic_polar_has_higher_moment(self): + # Alternating hydrophobic/polar gives high periodicity + # e.g. KALALALA at 100deg/residue should show amphipathic character + # vs uniform KKKKKKKKwhich is all charged + seq_amphipathic = "KLKLKLKL" + seq_uniform = "KKKKKKKK" + result_amph = hydrophobic_moment(seq_amphipathic) + result_unif = hydrophobic_moment(seq_uniform) + assert isinstance(result_amph, float) + assert isinstance(result_unif, float) + + def test_known_amp_has_nonzero_moment(self): + # KWKLFKKIGAVLKVL is a classic AMP (magainin analogue) + result = hydrophobic_moment("KWKLFKKIGAVLKVL") + assert result > 0.0 + + def test_returns_float_rounded_to_4dp(self): + result = hydrophobic_moment("KWKLFKK") + assert isinstance(result, float) + # Check 4 decimal places + assert result == round(result, 4) + + def test_single_residue(self): + result = hydrophobic_moment("K") + # sin(0) = 0, cos(0) = 1 → moment = |H_K * cos(0)| / 1 = |H_K| + assert isinstance(result, float) + + +class TestComputeFeaturesAmphipathicity: + def test_hydrophobic_moment_in_features(self): + features = compute_features("KWKLFKKIGAVLKVL") + assert "hydrophobic_moment" in features + assert isinstance(features["hydrophobic_moment"], float) + assert features["hydrophobic_moment"] >= 0.0 + + def test_empty_sequence_does_not_crash(self): + features = compute_features("") + assert "hydrophobic_moment" in features + assert features["hydrophobic_moment"] == 0.0 + + def test_all_canonical_amino_acids(self): + seq = "ACDEFGHIKLMNPQRSTVWY" + features = compute_features(seq) + assert "hydrophobic_moment" in features + assert features["hydrophobic_moment"] >= 0.0 diff --git a/tests/test_pipeline_filters.py b/tests/test_pipeline_filters.py new file mode 100644 index 00000000..a078063c --- /dev/null +++ b/tests/test_pipeline_filters.py @@ -0,0 +1,131 @@ +"""Tests for pipeline length filters and selection thresholds.""" +from __future__ import annotations + +import json +from pathlib import Path + +import pytest + +from openamp_foundry.pipeline import _passes_length_filter, run_ranking_pipeline, score_candidates + + +def test_length_filter_within_range(): + assert _passes_length_filter("KWKLFKKIGAVLKVL", min_length=8, max_length=35) is True + + +def test_length_filter_too_short(): + assert _passes_length_filter("KWKL", min_length=8, max_length=35) is False + + +def test_length_filter_too_long(): + assert _passes_length_filter("K" * 40, min_length=8, max_length=35) is False + + +def test_length_filter_boundary_min(): + assert _passes_length_filter("K" * 8, min_length=8, max_length=35) is True + + +def test_length_filter_boundary_max(): + assert _passes_length_filter("K" * 35, min_length=8, max_length=35) is True + + +def test_invalid_length_candidate_has_zero_activity(tmp_path): + """A candidate outside length bounds should score 0 for activity.""" + csv = tmp_path / "candidates.csv" + csv.write_text("id,sequence,source\nTEST-001,CCCCCCCCCCCC,test\n") + # CCCCCCCCCCCC is 12 chars but all Cys — valid alphabet but safety should be poor + scored, config = score_candidates(csv, config_path="configs/pipeline.yaml") + assert len(scored) == 1 + # All-cys peptide has non-zero safety risk + assert scored[0].scores["safety"] < 1.0 + + +def test_all_repeat_candidate_is_valid_but_has_low_charge(tmp_path): + """AAAAAAAAGGGGGGGG is valid amino acids but lacks positive charge (bad for AMPs).""" + csv = tmp_path / "candidates.csv" + csv.write_text("id,sequence,source\nTEST-001,AAAAAAAAGGGGGGGG,test\n") + scored, config = score_candidates(csv, config_path="configs/pipeline.yaml") + assert len(scored) == 1 + item = scored[0] + assert item.valid is True + # Net charge = 0, so charge_density = 0 → low charge score → low ensemble + assert item.features["net_charge_proxy"] == 0 + # Should score lower than a genuine AMP-like sequence (KWK-type) + assert item.scores["activity"] < 0.6 + + +def test_pipeline_marks_invalid_sequence(tmp_path): + """A sequence with non-canonical amino acids (B, X) should be marked invalid.""" + csv = tmp_path / "candidates.csv" + csv.write_text("id,sequence,source\nBAD-001,KWBXLFKK,test\n") + scored, _ = score_candidates(csv, config_path="configs/pipeline.yaml") + assert len(scored) == 1 + assert scored[0].valid is False + assert scored[0].scores["activity"] == 0.0 + assert any("non-canonical" in f for f in scored[0].known_failure_modes) + + +def test_selected_field_in_jsonl_output(tmp_path): + """The JSONL output should include a 'selected' boolean per candidate.""" + out = tmp_path / "ranked.jsonl" + run_ranking_pipeline( + candidate_path="examples/sequences/demo_candidates.csv", + reference_path="examples/known_reference/demo_known_amps.csv", + out_path=out, + ) + rows = [json.loads(line) for line in out.read_text().splitlines() if line.strip()] + assert all("selected" in row for row in rows) + assert any(row["selected"] for row in rows) + + +def test_run_manifest_generated(tmp_path): + """Run manifest should be auto-generated alongside the output.""" + out = tmp_path / "ranked.jsonl" + run_ranking_pipeline( + candidate_path="examples/sequences/demo_candidates.csv", + reference_path="examples/known_reference/demo_known_amps.csv", + out_path=out, + ) + manifest = tmp_path / "run_manifest.json" + assert manifest.exists() + data = json.loads(manifest.read_text()) + assert "run_id" in data + assert "pipeline_version" in data + assert "config_hash" in data + assert "input_hashes" in data + assert len(data["inputs"]) >= 1 + + +def test_run_manifest_explicit_path(tmp_path): + """Run manifest can be written to an explicit path.""" + out = tmp_path / "ranked.jsonl" + manifest_out = tmp_path / "my_manifest.json" + run_ranking_pipeline( + candidate_path="examples/sequences/demo_candidates.csv", + reference_path="examples/known_reference/demo_known_amps.csv", + out_path=out, + manifest_path=manifest_out, + ) + assert manifest_out.exists() + + +def test_novelty_filter_excludes_near_duplicates(tmp_path): + """Candidates with novelty below min_novelty threshold should not be selected.""" + csv = tmp_path / "candidates.csv" + refs = tmp_path / "refs.csv" + # Identical to reference — novelty will be 0.0 + csv.write_text("id,sequence,source\nTEST-001,KWKLFKKIGAVLKVL,test\n") + refs.write_text("id,sequence,source\nREF-001,KWKLFKKIGAVLKVL,reference\n") + out = tmp_path / "ranked.jsonl" + run_ranking_pipeline( + candidate_path=csv, + reference_path=refs, + out_path=out, + cert_dir=tmp_path / "certs", + ) + rows = [json.loads(line) for line in out.read_text().splitlines() if line.strip()] + # With novelty=0.0 and min_novelty=0.20 from config, should not be selected + assert not any(row["selected"] for row in rows) + # But should still be scored and appear in output + assert len(rows) == 1 + assert rows[0]["scores"]["novelty"] == 0.0