diff --git a/Makefile b/Makefile index 3a4679e0..11358f51 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -.PHONY: demo test lint clean bench-leakage bench-baseline bench-hidden-active generate phase3 pilot validate-scoring validate-scoring-strict external-predict pilot-confident presynth-qc gold-standard +.PHONY: demo test lint clean bench-leakage bench-baseline bench-hidden-active generate phase3 pilot validate-scoring validate-scoring-strict external-predict pilot-confident presynth-qc gold-standard diversity PYTHON := $(shell [ -f .venv/bin/python ] && echo .venv/bin/python || echo python3) PYTEST := $(shell [ -f .venv/bin/pytest ] && echo .venv/bin/pytest || echo pytest) @@ -109,5 +109,10 @@ gold-standard: --out outputs/gold_standard_calibration.md \ --config configs/pipeline.yaml +diversity: + PYTHONPATH=src $(PYTHON) -m openamp_foundry.cli diversity-check \ + --panel-csv outputs/confident_panel.csv \ + --out outputs/diversity_report.md + clean: rm -rf outputs/*.jsonl outputs/*.md outputs/*.json outputs/evidence outputs/phase3_evidence .pytest_cache .ruff_cache diff --git a/src/openamp_foundry/analysis/__init__.py b/src/openamp_foundry/analysis/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/src/openamp_foundry/analysis/diversity.py b/src/openamp_foundry/analysis/diversity.py new file mode 100644 index 00000000..352de891 --- /dev/null +++ b/src/openamp_foundry/analysis/diversity.py @@ -0,0 +1,183 @@ +"""Sequence diversity analysis for AMP candidate panels. + +Two concerns before synthesis: +1. Redundancy — candidates >60% similar to another in the panel occupy the same + structural space; testing both adds cost without proportional information gain. +2. Family-level risk — if all variants of a scaffold share a structural liability + (e.g., 5 trypsin sites in all SEED-005 14-mers), the whole family may fail + together, wiping out the corresponding synthesis investment. + +Key outputs: +- Pairwise similarity matrix (Levenshtein-based, no alignment required) +- Cluster assignments (single-linkage at configurable threshold) +- Recommended minimal diverse panel: best 1 candidate per cluster +- Family-level structural warnings +""" +from __future__ import annotations + + +def levenshtein_similarity(a: str, b: str) -> float: + """Sequence similarity in [0, 1] via normalised Levenshtein distance. + + Defined as 1 - edit_distance / max(len(a), len(b)). + Pure Python, O(|a| * |b|) — fine for short peptides. + """ + a, b = a.upper(), b.upper() + la, lb = len(a), len(b) + if la == 0 and lb == 0: + return 1.0 + dp = list(range(lb + 1)) + for ca in a: + ndp = [dp[0] + 1] + for j, cb in enumerate(b): + ndp.append(min(dp[j + 1] + 1, ndp[j] + 1, dp[j] + (0 if ca == cb else 1))) + dp = ndp + return 1.0 - dp[lb] / max(la, lb) + + +def pairwise_similarity_matrix( + sequences: list[str], +) -> list[list[float]]: + """Return an n×n symmetric similarity matrix.""" + n = len(sequences) + mat = [[0.0] * n for _ in range(n)] + for i in range(n): + mat[i][i] = 1.0 + for j in range(i + 1, n): + sim = levenshtein_similarity(sequences[i], sequences[j]) + mat[i][j] = sim + mat[j][i] = sim + return mat + + +def cluster_panel( + candidates: list[dict], + similarity_threshold: float = 0.60, +) -> list[dict]: + """Assign cluster IDs using greedy single-linkage clustering. + + Each candidate dict must have 'candidate_id' and 'sequence'. + Returns candidates with an added 'cluster_id' field. + Earlier (higher-priority) candidates seed clusters. + """ + sequences = [c["sequence"] for c in candidates] + cluster_ids: list[int] = [-1] * len(candidates) + next_cluster = 0 + + for i, seq_i in enumerate(sequences): + if cluster_ids[i] != -1: + continue + cluster_ids[i] = next_cluster + for j in range(i + 1, len(sequences)): + if cluster_ids[j] != -1: + continue + if levenshtein_similarity(seq_i, sequences[j]) >= similarity_threshold: + cluster_ids[j] = next_cluster + next_cluster += 1 + + return [ + {**c, "cluster_id": cluster_ids[i]} + for i, c in enumerate(candidates) + ] + + +def recommend_minimal_diverse_panel( + clustered: list[dict], + n_per_cluster: int = 1, +) -> list[dict]: + """Pick the top `n_per_cluster` candidates per cluster by list order. + + Assumes the input is already sorted by priority (highest first). + Returns the minimal set that covers all structural families. + """ + seen: dict[int, int] = {} + minimal: list[dict] = [] + for c in clustered: + cid = c["cluster_id"] + count = seen.get(cid, 0) + if count < n_per_cluster: + minimal.append(c) + seen[cid] = count + 1 + return minimal + + +def diversity_stats(clustered: list[dict]) -> dict: + """Compute diversity statistics for a clustered panel.""" + from collections import Counter + cluster_sizes = Counter(c["cluster_id"] for c in clustered) + n_clusters = len(cluster_sizes) + n_candidates = len(clustered) + redundant = sum(1 for c in clustered if cluster_sizes[c["cluster_id"]] > 1) + singletons = sum(1 for s in cluster_sizes.values() if s == 1) + largest_cluster = max(cluster_sizes.values()) + sequences = [c["sequence"] for c in clustered] + mat = pairwise_similarity_matrix(sequences) + # Mean off-diagonal similarity + n = len(sequences) + if n > 1: + total = sum(mat[i][j] for i in range(n) for j in range(n) if i != j) + mean_sim = total / (n * (n - 1)) + else: + mean_sim = 1.0 + return { + "n_candidates": n_candidates, + "n_clusters": n_clusters, + "n_redundant": redundant, + "n_singletons": singletons, + "largest_cluster_size": largest_cluster, + "mean_pairwise_similarity": round(mean_sim, 4), + "diversity_score": round(1.0 - mean_sim, 4), + } + + +def family_structural_warnings( + candidates: list[dict], + min_family_size: int = 3, +) -> list[dict]: + """Detect families where ALL members share a structural liability. + + Expects each candidate dict to optionally include: + - 'seed': e.g. 'SEED-005' + - 'trypsin_sites' (list) from presynth QC + - 'mu_h' (float) + - 'methionine_count' (int) + + Returns a list of warning dicts. + """ + from collections import defaultdict + families: dict[str, list[dict]] = defaultdict(list) + for c in candidates: + seed = c.get("seed", "UNKNOWN") + families[seed].append(c) + + warnings = [] + for seed, members in families.items(): + if len(members) < min_family_size: + continue + # Trypsin sites warning: all members have >= 4 sites + trypsin_counts = [len(c.get("trypsin_sites", [])) for c in members if "trypsin_sites" in c] + if trypsin_counts and all(t >= 4 for t in trypsin_counts): + warnings.append({ + "family": seed, + "n_members": len(members), + "warning_type": "TRYPSIN_STABILITY", + "message": ( + f"All {len(members)} {seed} members have ≥4 trypsin sites — " + "serum stability < 2h expected for the ENTIRE family. " + "If goal is in-vivo use, consider N-methylation or D-amino acid substitutions." + ), + }) + # High μH warning: all members μH > 0.65 + mu_h_vals = [c.get("mu_h", 0.0) for c in members if "mu_h" in c] + if mu_h_vals and all(m > 0.65 for m in mu_h_vals): + warnings.append({ + "family": seed, + "n_members": len(members), + "warning_type": "HEMOLYTIC_FAMILY_RISK", + "message": ( + f"All {len(members)} {seed} members have μH > 0.65 — " + "entire family is in the moderate-to-high hemolytic risk tier. " + "Test hemolysis before MIC; if HC50 < 2×MIC, family may lack a therapeutic window." + ), + }) + return warnings diff --git a/src/openamp_foundry/benchmark/evaluate.py b/src/openamp_foundry/benchmark/evaluate.py index 0f150352..e00ffbbf 100644 --- a/src/openamp_foundry/benchmark/evaluate.py +++ b/src/openamp_foundry/benchmark/evaluate.py @@ -29,14 +29,16 @@ def recall_at_k( def random_recall_at_k(n_candidates: int, n_positives: int, k: int) -> float: - """Expected recall@k for a random ranker. + """Expected recall@k for a random ranker (hypergeometric expectation). - E[recall@k] = min(k, n_positives) / n_positives when sampling without replacement. + E[hits in top-k] = k * n_positives / n_candidates + E[recall@k] = E[hits] / n_positives = k / n_candidates + + Capped at 1.0 when k >= n_candidates (all positives recovered). """ 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) + return round(min(1.0, k / n_candidates), 4) def enrichment_factor( diff --git a/src/openamp_foundry/cli.py b/src/openamp_foundry/cli.py index d91bd39d..fe8417e3 100644 --- a/src/openamp_foundry/cli.py +++ b/src/openamp_foundry/cli.py @@ -249,6 +249,37 @@ def build_parser() -> argparse.ArgumentParser: help="Pipeline config (for scoring weights).", ) + diversity_check = sub.add_parser( + "diversity-check", + help=( + "Analyse sequence redundancy in the confident panel. " + "Clusters candidates by sequence similarity and recommends a minimal " + "diverse synthesis set to maximise structural coverage per synthesis dollar." + ), + ) + diversity_check.add_argument( + "--panel-csv", + required=True, + help="Confident panel CSV (output of 'pilot-confident' command).", + ) + diversity_check.add_argument( + "--similarity-threshold", + type=float, + default=0.60, + help="Levenshtein similarity threshold for clustering (default: 0.60).", + ) + diversity_check.add_argument( + "--n-per-cluster", + type=int, + default=1, + help="Max candidates to keep per cluster in the minimal panel (default: 1).", + ) + diversity_check.add_argument( + "--out", + default="outputs/diversity_report.md", + help="Output markdown report path.", + ) + batch_pack = sub.add_parser( "batch-pack", help=( @@ -331,6 +362,9 @@ def main(argv: list[str] | None = None) -> int: if args.command == "gold-standard": return _run_gold_standard(args) + if args.command == "diversity-check": + return _run_diversity_check(args) + parser.error("unknown command") return 2 @@ -847,5 +881,256 @@ def _score(seq: str) -> dict: return 0 +def _run_diversity_check(args: argparse.Namespace) -> int: + import csv as _csv + from openamp_foundry.analysis.diversity import ( + cluster_panel, + diversity_stats, + family_structural_warnings, + pairwise_similarity_matrix, + recommend_minimal_diverse_panel, + ) + from openamp_foundry.features.physchem import compute_features + from openamp_foundry.qc.presynth_check import check_panel + + # Load panel + with open(args.panel_csv, newline="", encoding="utf-8") as f: + panel = list(_csv.DictReader(f)) + + # Pre-synthesis QC to get μH and trypsin sites + mu_h_map = {} + for row in panel: + seq = row["sequence"] + feats = compute_features(seq) + mu_h_map[row["candidate_id"]] = feats.get("hydrophobic_moment", 0.0) + + qc_results = check_panel(panel, mu_h_map=mu_h_map) + qc_by_id = {q.candidate_id: q for q in qc_results} + + # Cluster + threshold = args.similarity_threshold + clustered = cluster_panel(panel, similarity_threshold=threshold) + stats = diversity_stats(clustered) + minimal = recommend_minimal_diverse_panel(clustered, n_per_cluster=args.n_per_cluster) + + # Family-level warnings (merge QC data) + enriched = [] + for c in clustered: + qc = qc_by_id.get(c["candidate_id"]) + if qc: + enriched.append({ + **c, + "mu_h": qc.mu_h, + "trypsin_sites": qc.trypsin_sites, + "methionine_count": qc.methionine_count, + }) + else: + enriched.append(c) + fam_warnings = family_structural_warnings(enriched, min_family_size=3) + + # Pairwise matrix for report + sequences = [c["sequence"] for c in panel] + ids = [c["candidate_id"] for c in panel] + mat = pairwise_similarity_matrix(sequences) + + # Build report + from datetime import datetime, timezone + ts = datetime.now(timezone.utc).isoformat() + + lines = [ + "# Sequence Diversity Analysis", + "", + f"> Generated: {ts} ", + f"> Panel: {args.panel_csv} ", + f"> Similarity threshold: {threshold} ", + "", + "## Summary", + "", + "| Metric | Value |", + "|---|---|", + f"| Candidates analysed | {stats['n_candidates']} |", + f"| Structural clusters | {stats['n_clusters']} |", + f"| Redundant candidates (same cluster as earlier) | {stats['n_redundant']} |", + f"| Singleton clusters (fully unique) | {stats['n_singletons']} |", + f"| Largest cluster size | {stats['largest_cluster_size']} |", + f"| Mean pairwise similarity | {stats['mean_pairwise_similarity']:.3f} |", + f"| Diversity score (1 − mean_sim) | {stats['diversity_score']:.3f} |", + "", + ] + + if stats["n_redundant"] > 0: + lines += [ + "## ⚠ Redundancy Warning", + "", + f"{stats['n_redundant']} of {stats['n_candidates']} candidates share a cluster " + f"with an earlier (higher-priority) candidate at similarity threshold {threshold}.", + "Synthesising all of them adds cost with diminishing structural information gain.", + "", + ] + + # Cluster assignments table + lines += [ + "## Cluster Assignments", + "", + "| Rank | ID | Sequence | Cluster | Ensemble | Safety | μH | Status |", + "|--:|---|---|:---:|:---:|:---:|:---:|:---:|", + ] + + minimal_ids = {c["candidate_id"] for c in minimal} + seen_clusters: set[int] = set() + for c in clustered: + cid = c["cluster_id"] + seen_clusters.add(cid) + rank = c.get("pilot_rank", "") + cand_id = c["candidate_id"] + seq = c["sequence"] + ens = float(c.get("ensemble", 0)) + safe = float(c.get("safety", 0)) + mu_h = mu_h_map.get(cand_id, 0.0) + status = "✅ KEEP" if cand_id in minimal_ids else "⚠ REDUNDANT" + lines.append( + f"| {rank} | {cand_id} | `{seq}` | C{cid} | {ens:.3f} | {safe:.3f} | {mu_h:.2f} | {status} |" + ) + + # Minimal diverse panel + lines += [ + "", + f"## Recommended Minimal Diverse Panel ({len(minimal)} candidates)", + "", + f"Best 1 candidate per structural cluster. Covers {stats['n_clusters']} distinct " + f"structural families in {len(minimal)} synthesis slots instead of {stats['n_candidates']}.", + "", + "| Priority | ID | Sequence | Cluster | Ensemble | Safety | μH |", + "|:---:|---|---|:---:|:---:|:---:|:---:|", + ] + for i, c in enumerate(minimal, 1): + cand_id = c["candidate_id"] + seq = c["sequence"] + ens = float(c.get("ensemble", 0)) + safe = float(c.get("safety", 0)) + mu_h = mu_h_map.get(cand_id, 0.0) + cluster_id = c["cluster_id"] + lines.append( + f"| {i} | {cand_id} | `{seq}` | C{cluster_id} | {ens:.3f} | {safe:.3f} | {mu_h:.2f} |" + ) + + # Optimal cluster representatives (best ensemble per cluster, tie-break by low μH + no Met) + from collections import defaultdict as _dd + by_cluster: dict = _dd(list) + for c in clustered: + by_cluster[c["cluster_id"]].append(c) + + optimal = [] + for cid in sorted(by_cluster.keys()): + members = by_cluster[cid] + # Sort: highest ensemble first, tie-break: no Met then lowest μH + members_sorted = sorted( + members, + key=lambda x: ( + -float(x.get("ensemble", 0)), + 1 if qc_by_id.get(x["candidate_id"]) and qc_by_id[x["candidate_id"]].methionine_count > 0 else 0, + mu_h_map.get(x["candidate_id"], 0.0), + ) + ) + for rep in members_sorted[:args.n_per_cluster]: + optimal.append(rep) + + optimal_ids = {c["candidate_id"] for c in optimal} + if optimal_ids != minimal_ids: + lines += [ + "", + f"## Optimal Cluster Representatives ({len(optimal)} candidates)", + "", + "Picked by **highest ensemble score per cluster** (tie-break: no Met, lower μH). " + "This differs from 'Recommended Minimal' above, which picked by pilot rank.", + "", + "| Priority | ID | Sequence | Cluster | Ensemble | Safety | μH | Met |", + "|:---:|---|---|:---:|:---:|:---:|:---:|:---:|", + ] + for i, c in enumerate(optimal, 1): + cand_id = c["candidate_id"] + seq = c["sequence"] + ens = float(c.get("ensemble", 0)) + safe = float(c.get("safety", 0)) + mu_h = mu_h_map.get(cand_id, 0.0) + cluster_id = c["cluster_id"] + qc = qc_by_id.get(cand_id) + met = qc.methionine_count if qc else 0 + met_flag = f"⚠ {met}×Met" if met > 0 else "✓" + lines.append( + f"| {i} | {cand_id} | `{seq}` | C{cluster_id} | {ens:.3f} | {safe:.3f} | {mu_h:.2f} | {met_flag} |" + ) + lines += [ + "", + "> Use these as the actual synthesis candidates if ordering ≤4 peptides.", + ] + + # Family-level warnings + if fam_warnings: + lines += ["", "## Family-Level Structural Warnings", ""] + for w in fam_warnings: + lines += [f"**{w['family']} — {w['warning_type']}:** {w['message']}", ""] + + # Pairwise similarity heatmap (abbreviated) + lines += [ + "", + "## Pairwise Similarity (upper triangle, threshold shown)", + "", + "Values ≥ threshold are flagged as same-cluster.", + "", + "| | " + " | ".join(f"**{i}**" for i in range(1, len(ids) + 1)) + " |", + "|---|" + "---|" * len(ids), + ] + for i, row_id in enumerate(ids): + cells = [] + for j in range(len(ids)): + if j < i: + cells.append("") + elif j == i: + cells.append("1.00") + else: + v = mat[i][j] + cells.append(f"**{v:.2f}**" if v >= threshold else f"{v:.2f}") + lines.append(f"| **{i+1}** | " + " | ".join(cells) + " |") + + lines += [ + "", + "## Synthesis Recommendation", + "", + f"**If budget allows ≤{len(minimal)} peptides:** Use the minimal diverse panel above.", + f"**If budget allows all {stats['n_candidates']} peptides:** Order all — redundancy " + "provides SAR (structure-activity relationship) data within each scaffold family.", + "", + "The minimal panel is strongly preferred for a first-pass experiment:", + "- Covers all structural families", + "- Maximises the chance of finding at least one hit in each scaffold class", + "- Preserves budget for a second wave if any scaffold family shows promise", + "", + "## Disclaimer", + "", + "Sequence similarity (Levenshtein) is a proxy for structural similarity.", + "Peptides with different sequences can still adopt similar 3D conformations,", + "and vice versa. For short AMPs (<20 AA), sequence similarity is generally a", + "reliable proxy for functional similarity.", + ] + + out = Path(args.out) + out.parent.mkdir(parents=True, exist_ok=True) + out.write_text("\n".join(lines) + "\n", encoding="utf-8") + + print(json.dumps({ + "status": "ok", + "n_candidates": stats["n_candidates"], + "n_clusters": stats["n_clusters"], + "n_redundant": stats["n_redundant"], + "minimal_panel_size": len(minimal), + "mean_pairwise_similarity": stats["mean_pairwise_similarity"], + "diversity_score": stats["diversity_score"], + "n_family_warnings": len(fam_warnings), + "out": args.out, + }, indent=2)) + return 0 + + if __name__ == "__main__": raise SystemExit(main()) diff --git a/src/openamp_foundry/qc/presynth_check.py b/src/openamp_foundry/qc/presynth_check.py index bc227ab0..bf9d4681 100644 --- a/src/openamp_foundry/qc/presynth_check.py +++ b/src/openamp_foundry/qc/presynth_check.py @@ -44,8 +44,6 @@ def _net_charge_at_ph(seq: str, ph: float) -> float: charge -= 1.0 / (1.0 + 10 ** (_PKA["C_term"] - ph)) # Basic residues (positive) for aa, pka in [("K", _PKA["K"]), ("R", _PKA["R"]), ("H", _PKA["H"])]: - for _ in seq: - pass charge += seq.count(aa) / (1.0 + 10 ** (ph - pka)) # Acidic residues (negative) for aa, pka in [("D", _PKA["D"]), ("E", _PKA["E"]), ("C", _PKA["C"]), ("Y", _PKA["Y"])]: @@ -166,7 +164,10 @@ def check_sequence(candidate_id: str, seq: str, mu_h: float = 0.0) -> SynthQC: qc.hydrophobic_run = max(runs, key=len) qc.aggregation_risk = True - # Proteolytic sites (interior only — terminal K/R are typical, less risky) + # Proteolytic sites (interior only — terminal K/R are typical, less risky). + # _TRYPSIN_RE already requires a following character via (?=.), so the + # m.start() < len(seq)-1 check is the canonical guard; the regex lookahead + # provides a defence-in-depth backstop. qc.trypsin_sites = [m.start() for m in _TRYPSIN_RE.finditer(seq) if m.start() < len(seq) - 1] qc.chymotrypsin_sites = [m.start() for m in _CHYMOTRYPSIN_RE.finditer(seq) if m.start() < len(seq) - 1] qc.deamidation_sites = [ @@ -239,11 +240,15 @@ def check_panel( ) -> list[SynthQC]: """Run QC on a list of candidate dicts (must have 'candidate_id' and 'sequence').""" mu_h_map = mu_h_map or {} - return [ - check_sequence( - c["candidate_id"], - c["sequence"], - mu_h=mu_h_map.get(c["candidate_id"], 0.0), - ) - for c in candidates - ] + results = [] + for c in candidates: + try: + cid = c["candidate_id"] + seq = c["sequence"] + except KeyError as exc: + raise KeyError( + f"Candidate row missing required column {exc} — " + "check that your panel CSV has 'candidate_id' and 'sequence' headers." + ) from exc + results.append(check_sequence(cid, seq, mu_h=mu_h_map.get(cid, 0.0))) + return results diff --git a/src/openamp_foundry/scoring/activity.py b/src/openamp_foundry/scoring/activity.py index 615e9215..52dc3c48 100644 --- a/src/openamp_foundry/scoring/activity.py +++ b/src/openamp_foundry/scoring/activity.py @@ -42,10 +42,13 @@ def activity_likeness_score(features: dict) -> float: # 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 + # Base weights (0.24+0.27+0.17 = 0.68) leave 0.32 headroom for bonuses + # (aromatic max 0.10 + amphipathicity max 0.15 = 0.25), ceiling = 0.93 < 1.0. + # Scale factor vs original (0.28, 0.32, 0.20): multiply each by 0.68/0.80 ≈ 0.85. score = ( - 0.28 * length_score - + 0.32 * charge_score - + 0.20 * hydro_score + 0.24 * length_score + + 0.27 * charge_score + + 0.17 * hydro_score + aromatic_bonus + amphipathicity_score ) diff --git a/tests/test_benchmark_evaluate.py b/tests/test_benchmark_evaluate.py index 3e6c7a5a..c6342c05 100644 --- a/tests/test_benchmark_evaluate.py +++ b/tests/test_benchmark_evaluate.py @@ -58,10 +58,20 @@ def test_empty_positives_returns_zero(self): 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 + # E[recall@k] = k / n_candidates (hypergeometric expectation) + # k=2, n_candidates=5 → 2/5 = 0.4 result = random_recall_at_k(n_candidates=5, n_positives=2, k=2) - assert 0 < result < 1 + assert abs(result - 0.4) < 0.001, f"Expected 0.4, got {result}" + + def test_expected_random_recall_k_less_than_positives(self): + # k=1 < n_positives=5, n_candidates=10 → E[recall@1] = 1/10 = 0.1 + result = random_recall_at_k(n_candidates=10, n_positives=5, k=1) + assert abs(result - 0.1) < 0.001, f"Expected 0.1, got {result}" + + def test_random_recall_real_benchmark_scale(self): + # 44 AMPs, 88 total candidates, k=8 → E[recall@8] = 8/88 ≈ 0.0909 + result = random_recall_at_k(n_candidates=88, n_positives=44, k=8) + assert abs(result - 8 / 88) < 0.001, f"Expected {8/88:.4f}, got {result}" def test_zero_candidates_returns_zero(self): assert random_recall_at_k(0, 2, 2) == 0.0 @@ -97,8 +107,7 @@ def test_summary_has_required_keys(self): 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() + assert "do not prove biological efficacy" 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 diff --git a/tests/test_diversity.py b/tests/test_diversity.py new file mode 100644 index 00000000..1d0f292c --- /dev/null +++ b/tests/test_diversity.py @@ -0,0 +1,291 @@ +"""Tests for src/openamp_foundry/analysis/diversity.py""" +import pytest +from openamp_foundry.analysis.diversity import ( + cluster_panel, + diversity_stats, + family_structural_warnings, + levenshtein_similarity, + pairwise_similarity_matrix, + recommend_minimal_diverse_panel, +) + + +# --------------------------------------------------------------------------- +# levenshtein_similarity +# --------------------------------------------------------------------------- + +class TestLevenshteinSimilarity: + def test_identical_sequences(self): + assert levenshtein_similarity("ACDEF", "ACDEF") == 1.0 + + def test_completely_different(self): + # Every character substituted + sim = levenshtein_similarity("AAAAA", "BBBBB") + assert sim == 0.0 + + def test_empty_both(self): + assert levenshtein_similarity("", "") == 1.0 + + def test_single_substitution(self): + # 1 edit in 5 → similarity = 1 - 1/5 = 0.8 + sim = levenshtein_similarity("ACDEF", "ACGEF") + assert abs(sim - 0.8) < 1e-9 + + def test_one_insertion(self): + # "ABC" → "ABCD": 1 insertion, max_len = 4 → sim = 1 - 1/4 = 0.75 + sim = levenshtein_similarity("ABC", "ABCD") + assert abs(sim - 0.75) < 1e-9 + + def test_symmetric(self): + a, b = "KRLFKKIG", "KRLFKKIGSALKFL" + assert levenshtein_similarity(a, b) == levenshtein_similarity(b, a) + + def test_case_insensitive(self): + assert levenshtein_similarity("acdef", "ACDEF") == 1.0 + + def test_known_seed005_variants_high_similarity(self): + # Two SEED-005 variants — should be > 0.6 + v1 = "KRLFKKIGSALKFL" + v2 = "KRLFKKVGSALRFL" + assert levenshtein_similarity(v1, v2) > 0.6 + + def test_different_families_low_similarity(self): + seed005 = "KRLFKKIGSALKFL" + seed003 = "RRWQWRMKKLG" + assert levenshtein_similarity(seed005, seed003) < 0.5 + + +# --------------------------------------------------------------------------- +# pairwise_similarity_matrix +# --------------------------------------------------------------------------- + +class TestPairwiseSimilarityMatrix: + def test_diagonal_is_one(self): + seqs = ["ACDEF", "KLMNO", "PQRST"] + mat = pairwise_similarity_matrix(seqs) + for i in range(len(seqs)): + assert mat[i][i] == 1.0 + + def test_symmetric(self): + seqs = ["ACDEF", "KLMNO", "PQRST"] + mat = pairwise_similarity_matrix(seqs) + n = len(seqs) + for i in range(n): + for j in range(n): + assert abs(mat[i][j] - mat[j][i]) < 1e-9 + + def test_shape(self): + seqs = ["ABC", "DEF", "GHI", "JKL"] + mat = pairwise_similarity_matrix(seqs) + assert len(mat) == 4 + assert all(len(row) == 4 for row in mat) + + def test_single_sequence(self): + mat = pairwise_similarity_matrix(["ACDEF"]) + assert mat == [[1.0]] + + def test_values_in_range(self): + seqs = ["KRLFKKIGSALKFL", "KRLFKKVGSALRFL", "RRWQWRMKKLG"] + mat = pairwise_similarity_matrix(seqs) + for row in mat: + for v in row: + assert 0.0 <= v <= 1.0 + + +# --------------------------------------------------------------------------- +# cluster_panel +# --------------------------------------------------------------------------- + +class TestClusterPanel: + def _make_candidates(self, id_seq_pairs): + return [{"candidate_id": cid, "sequence": seq} for cid, seq in id_seq_pairs] + + def test_identical_sequences_same_cluster(self): + candidates = self._make_candidates([ + ("A", "KRLFKKIGSALKFL"), + ("B", "KRLFKKIGSALKFL"), + ]) + clustered = cluster_panel(candidates, similarity_threshold=0.6) + assert clustered[0]["cluster_id"] == clustered[1]["cluster_id"] + + def test_different_families_different_clusters(self): + candidates = self._make_candidates([ + ("seed005", "KRLFKKIGSALKFL"), + ("seed003", "RRWQWRMKKLG"), + ]) + clustered = cluster_panel(candidates, similarity_threshold=0.6) + assert clustered[0]["cluster_id"] != clustered[1]["cluster_id"] + + def test_preserves_original_fields(self): + candidates = [{"candidate_id": "X", "sequence": "ACDEF", "ensemble": 0.85}] + clustered = cluster_panel(candidates) + assert clustered[0]["ensemble"] == 0.85 + assert "cluster_id" in clustered[0] + + def test_cluster_id_added(self): + candidates = self._make_candidates([("A", "ACDEF"), ("B", "KLMNO")]) + clustered = cluster_panel(candidates) + assert all("cluster_id" in c for c in clustered) + + def test_threshold_zero_all_same_cluster(self): + # threshold=0 → everything clusters together + candidates = self._make_candidates([ + ("A", "AAAA"), ("B", "BBBB"), ("C", "CCCC"), + ]) + clustered = cluster_panel(candidates, similarity_threshold=0.0) + ids = {c["cluster_id"] for c in clustered} + assert len(ids) == 1 + + def test_threshold_one_all_different(self): + # threshold=1.0 → only identical sequences cluster + candidates = self._make_candidates([ + ("A", "AAAA"), ("B", "BBBB"), ("C", "CCCC"), + ]) + clustered = cluster_panel(candidates, similarity_threshold=1.0) + ids = {c["cluster_id"] for c in clustered} + assert len(ids) == 3 + + def test_seed005_family_clusters_together(self): + candidates = self._make_candidates([ + ("V068", "KRLMKKIGSAIKFL"), + ("V009", "KRFFKKIGSALKFA"), + ("V063", "KRLFRKIGSALKFV"), + ]) + clustered = cluster_panel(candidates, similarity_threshold=0.6) + ids = {c["cluster_id"] for c in clustered} + # All similar → should share cluster + assert len(ids) == 1 + + +# --------------------------------------------------------------------------- +# recommend_minimal_diverse_panel +# --------------------------------------------------------------------------- + +class TestRecommendMinimalDiversePanel: + def test_one_per_cluster(self): + clustered = [ + {"candidate_id": "A", "cluster_id": 0}, + {"candidate_id": "B", "cluster_id": 0}, + {"candidate_id": "C", "cluster_id": 1}, + ] + minimal = recommend_minimal_diverse_panel(clustered, n_per_cluster=1) + assert len(minimal) == 2 + assert minimal[0]["candidate_id"] == "A" + assert minimal[1]["candidate_id"] == "C" + + def test_two_per_cluster(self): + clustered = [ + {"candidate_id": "A", "cluster_id": 0}, + {"candidate_id": "B", "cluster_id": 0}, + {"candidate_id": "C", "cluster_id": 0}, + {"candidate_id": "D", "cluster_id": 1}, + ] + minimal = recommend_minimal_diverse_panel(clustered, n_per_cluster=2) + ids = [c["candidate_id"] for c in minimal] + assert ids == ["A", "B", "D"] + + def test_all_unique_clusters(self): + clustered = [ + {"candidate_id": str(i), "cluster_id": i} for i in range(5) + ] + minimal = recommend_minimal_diverse_panel(clustered) + assert len(minimal) == 5 + + def test_preserves_fields(self): + clustered = [{"candidate_id": "A", "cluster_id": 0, "ensemble": 0.9}] + minimal = recommend_minimal_diverse_panel(clustered) + assert minimal[0]["ensemble"] == 0.9 + + +# --------------------------------------------------------------------------- +# diversity_stats +# --------------------------------------------------------------------------- + +class TestDiversityStats: + def test_all_unique(self): + clustered = [{"sequence": seq, "cluster_id": i} + for i, seq in enumerate(["ACDEF", "KLMNO", "PQRST"])] + stats = diversity_stats(clustered) + assert stats["n_clusters"] == 3 + assert stats["n_redundant"] == 0 + assert stats["n_singletons"] == 3 + + def test_all_same_cluster(self): + clustered = [{"sequence": "ACDEF", "cluster_id": 0} for _ in range(3)] + stats = diversity_stats(clustered) + assert stats["n_clusters"] == 1 + assert stats["n_redundant"] == 3 + assert stats["largest_cluster_size"] == 3 + + def test_diversity_score_range(self): + clustered = [ + {"sequence": "AAAA", "cluster_id": 0}, + {"sequence": "BBBB", "cluster_id": 1}, + ] + stats = diversity_stats(clustered) + assert 0.0 <= stats["diversity_score"] <= 1.0 + + def test_single_candidate(self): + clustered = [{"sequence": "ACDEF", "cluster_id": 0}] + stats = diversity_stats(clustered) + assert stats["n_candidates"] == 1 + assert stats["mean_pairwise_similarity"] == 1.0 + + +# --------------------------------------------------------------------------- +# family_structural_warnings +# --------------------------------------------------------------------------- + +class TestFamilyStructuralWarnings: + def _make_member(self, seed, trypsin_sites, mu_h, cid="X"): + return { + "candidate_id": cid, + "seed": seed, + "trypsin_sites": trypsin_sites, + "mu_h": mu_h, + } + + def test_trypsin_warning_fires_for_large_family(self): + members = [ + self._make_member("SEED-005", [0, 1, 4, 5, 11], 0.70, f"V{i}") + for i in range(3) + ] + warnings = family_structural_warnings(members, min_family_size=3) + types = [w["warning_type"] for w in warnings] + assert "TRYPSIN_STABILITY" in types + + def test_no_warning_for_small_family(self): + members = [ + self._make_member("SEED-005", [0, 1, 4], 0.70, f"V{i}") + for i in range(2) + ] + warnings = family_structural_warnings(members, min_family_size=3) + assert len(warnings) == 0 + + def test_hemolytic_warning_fires_when_all_high_mu_h(self): + members = [ + self._make_member("SEED-005", [], 0.75, f"V{i}") + for i in range(3) + ] + warnings = family_structural_warnings(members, min_family_size=3) + types = [w["warning_type"] for w in warnings] + assert "HEMOLYTIC_FAMILY_RISK" in types + + def test_no_hemolytic_warning_when_one_low_mu_h(self): + members = [ + self._make_member("SEED-005", [], 0.75, "V0"), + self._make_member("SEED-005", [], 0.80, "V1"), + self._make_member("SEED-005", [], 0.40, "V2"), # low μH + ] + warnings = family_structural_warnings(members, min_family_size=3) + types = [w["warning_type"] for w in warnings] + assert "HEMOLYTIC_FAMILY_RISK" not in types + + def test_warning_contains_family_name(self): + members = [ + self._make_member("SEED-005", [0, 1, 4, 5, 11], 0.80, f"V{i}") + for i in range(3) + ] + warnings = family_structural_warnings(members, min_family_size=3) + for w in warnings: + assert "SEED-005" in w["message"] diff --git a/tests/test_physchem_amphipathicity.py b/tests/test_physchem_amphipathicity.py index 0057f1ee..7439bba7 100644 --- a/tests/test_physchem_amphipathicity.py +++ b/tests/test_physchem_amphipathicity.py @@ -17,14 +17,16 @@ def test_uniform_sequence_low_moment(self): 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 + # e.g. KLKLKLKL at 100deg/residue should show amphipathic character + # vs KKKKKKKK which is all charged (low hydrophobic contrast) 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) + assert result_amph > result_unif, ( + f"Amphipathic KLKLKLKL moment ({result_amph:.4f}) should exceed " + f"uniform KKKKKKKK moment ({result_unif:.4f})" + ) def test_known_amp_has_nonzero_moment(self): # KWKLFKKIGAVLKVL is a classic AMP (magainin analogue) diff --git a/tests/test_presynth_qc.py b/tests/test_presynth_qc.py new file mode 100644 index 00000000..62420efc --- /dev/null +++ b/tests/test_presynth_qc.py @@ -0,0 +1,295 @@ +"""Tests for src/openamp_foundry/qc/presynth_check.py""" +import pytest +from openamp_foundry.qc.presynth_check import ( + SynthQC, + check_panel, + check_sequence, +) + + +# --------------------------------------------------------------------------- +# Molecular weight +# --------------------------------------------------------------------------- + +class TestMolecularWeight: + def test_glycine_monomer(self): + # Single G: residue mass 75.03 - 0 water corrections = 75.03 + qc = check_sequence("G1", "G") + assert qc.molecular_weight_da == pytest.approx(75.03, abs=0.5) + + def test_short_peptide_positive(self): + qc = check_sequence("test", "ACDEF") + assert qc.molecular_weight_da > 0 + + def test_longer_peptide_heavier(self): + qc_short = check_sequence("s", "KRLF") + qc_long = check_sequence("l", "KRLFKKIGSALKFL") + assert qc_long.molecular_weight_da > qc_short.molecular_weight_da + + def test_known_seed003(self): + # RRWQWRMKKLG should be ~1.4–1.6 kDa + qc = check_sequence("s3", "RRWQWRMKKLG") + assert 1300 < qc.molecular_weight_da < 1700 + + +# --------------------------------------------------------------------------- +# Isoelectric point +# --------------------------------------------------------------------------- + +class TestIsoelectricPoint: + def test_cationic_peptide_basic_pi(self): + # KKKK should be highly basic + qc = check_sequence("k4", "KKKK") + assert qc.isoelectric_point > 9.0 + + def test_acidic_peptide_low_pi(self): + # DDDD should be acidic + qc = check_sequence("d4", "DDDD") + assert qc.isoelectric_point < 5.0 + + def test_pi_in_valid_range(self): + qc = check_sequence("s", "KRLFKKIGSALKFL") + assert 0.0 < qc.isoelectric_point < 14.0 + + +# --------------------------------------------------------------------------- +# Net charge +# --------------------------------------------------------------------------- + +class TestNetCharge: + def test_cationic_amp_positive_charge(self): + # SEED-005 scaffold — K×3, R×1 → strongly positive + qc = check_sequence("s", "KRLFKKIGSALKFL") + assert qc.charge_ph74 > 2.0 + + def test_charge_decreases_at_higher_ph(self): + qc = check_sequence("s", "KRLFKKIGSALKFL") + assert qc.charge_ph74 < qc.charge_ph60 + + def test_lysine_rich_highly_positive(self): + qc = check_sequence("k", "KKKKKK") + assert qc.charge_ph74 > 4.0 + + def test_acidic_peptide_negative(self): + qc = check_sequence("d", "DDDDDD") + assert qc.charge_ph74 < 0.0 + + +# --------------------------------------------------------------------------- +# Oxidation risk (Met, Cys) +# --------------------------------------------------------------------------- + +class TestOxidationRisk: + def test_no_met_no_cys_no_risk(self): + qc = check_sequence("s", "KRLFKKIGSALKFL") + assert not qc.has_oxidation_risk + assert qc.methionine_count == 0 + assert qc.cysteine_count == 0 + + def test_met_triggers_oxidation_risk(self): + qc = check_sequence("s", "KRLMKKIGSAIKFL") + assert qc.has_oxidation_risk + assert qc.methionine_count == 1 + + def test_cys_triggers_oxidation_risk(self): + qc = check_sequence("s", "ACYCRIPACIAGER") + assert qc.has_oxidation_risk + assert qc.cysteine_count >= 2 + + def test_two_met_counted(self): + qc = check_sequence("s", "RRWNWRMKKMG") + assert qc.methionine_count == 2 + + def test_met_flag_in_flags_list(self): + qc = check_sequence("s", "KRLMKKIGSAIKFL") + met_flags = [f for f in qc.flags if "MET" in f] + assert len(met_flags) == 1 + + +# --------------------------------------------------------------------------- +# Trypsin/chymotrypsin sites +# --------------------------------------------------------------------------- + +class TestProteolytic: + def test_trypsin_sites_detected(self): + # KRLFKKIGSALKFL: K and R (not terminal, not before P) + qc = check_sequence("s", "KRLFKKIGSALKFL") + assert len(qc.trypsin_sites) >= 3 + + def test_trypsin_before_proline_not_counted(self): + # KP: K before P should NOT be counted + qc = check_sequence("s", "AKPGG") + # K at pos 1 (0-based) is before P → should NOT be in sites + # A at pos 0 is just a residue + assert 1 not in qc.trypsin_sites + + def test_chymotrypsin_w_detected(self): + qc = check_sequence("s", "RRWQWRMKKLG") + assert len(qc.chymotrypsin_sites) >= 1 + + def test_terminal_k_not_trypsin_site(self): + # K at last position should not be counted (terminal) + qc = check_sequence("s", "GGGGGK") + assert (len(qc.sequence) - 1) not in qc.trypsin_sites + + +# --------------------------------------------------------------------------- +# Deamidation +# --------------------------------------------------------------------------- + +class TestDeamidation: + def test_ng_hotspot_detected(self): + qc = check_sequence("s", "ANGGG") + assert any("N" in s for s in qc.deamidation_sites) + + def test_ns_hotspot_detected(self): + qc = check_sequence("s", "ANSGGG") + assert any("N" in s for s in qc.deamidation_sites) + + def test_no_hotspot_when_absent(self): + qc = check_sequence("s", "KRLFKKIGSALKFL") + assert qc.deamidation_sites == [] + + def test_hotspot_format(self): + qc = check_sequence("s", "ANGGG") + # Should be like "N1G" (residue + 1-based position + next residue) + assert len(qc.deamidation_sites) >= 1 + for site in qc.deamidation_sites: + assert site[0] == "N" + + +# --------------------------------------------------------------------------- +# Aggregation / hydrophobic runs +# --------------------------------------------------------------------------- + +class TestAggregationRisk: + def test_four_hydrophobic_run_flagged(self): + qc = check_sequence("s", "VILLKK") + assert qc.aggregation_risk + assert "VILL" in qc.hydrophobic_run + + def test_three_hydrophobic_not_flagged(self): + qc = check_sequence("s", "VILKK") + assert not qc.aggregation_risk + + def test_seed005_flagged_for_alkfl_run(self): + # KRLFKKIGSALKFL has LKFL at the end with ALKFL — let's check + qc = check_sequence("s", "KRLFKKIGSALKFL") + # "ALKFL" has 5 chars with A, L, K, F, L — K is not in VILMFW + # so not a pure hydrophobic run. Let's just check the function works + assert isinstance(qc.aggregation_risk, bool) + + def test_run_length_correct(self): + qc = check_sequence("s", "KVILMFWAAAK") + assert len(qc.hydrophobic_run) >= 4 + + +# --------------------------------------------------------------------------- +# UV chromophore +# --------------------------------------------------------------------------- + +class TestUVChromophore: + def test_w_gives_chromophore(self): + qc = check_sequence("s", "RRWQWRMKKLG") + assert qc.has_uv_chromophore + + def test_y_gives_chromophore(self): + qc = check_sequence("s", "KRWQYRMKKLG") + assert qc.has_uv_chromophore + + def test_f_gives_chromophore(self): + qc = check_sequence("s", "KRLFKKIGSALKFL") + assert qc.has_uv_chromophore + + def test_no_chromophore(self): + qc = check_sequence("s", "KRKRKLLL") + assert not qc.has_uv_chromophore + + def test_w_formulation_note_includes_eps280(self): + qc = check_sequence("s", "RRWWRR") + assert "A280" in qc.formulation_note + assert "ε=" in qc.formulation_note + + +# --------------------------------------------------------------------------- +# Hemolysis stratification +# --------------------------------------------------------------------------- + +class TestHemolysisStratification: + def test_low_mu_h_start_at_mic(self): + qc = check_sequence("s", "AAAAAAA", mu_h=0.40) + assert "Start at MIC" in qc.hemolysis_start_conc + + def test_moderate_mu_h_start_at_mic_third(self): + qc = check_sequence("s", "AAAAAAA", mu_h=0.65) + assert "MIC/3" in qc.hemolysis_start_conc + + def test_high_mu_h_start_at_mic_tenth(self): + qc = check_sequence("s", "AAAAAAA", mu_h=0.90) + assert "MIC/10" in qc.hemolysis_start_conc + + def test_boundary_at_055(self): + # μH exactly 0.55 → should still be "low risk" (> 0.55 triggers moderate) + qc = check_sequence("s", "AAAAAAA", mu_h=0.55) + assert "Start at MIC;" in qc.hemolysis_start_conc + + def test_boundary_at_080(self): + # μH exactly 0.80 → moderate (> 0.80 triggers high) + qc = check_sequence("s", "AAAAAAA", mu_h=0.80) + assert "MIC/3" in qc.hemolysis_start_conc + + +# --------------------------------------------------------------------------- +# Synthesis difficulty +# --------------------------------------------------------------------------- + +class TestSynthesisDifficulty: + def test_clean_peptide_low_difficulty(self): + # Simple 8-mer with no risky residues and no trypsin sites beyond terminal + qc = check_sequence("s", "RRWQWRMK") # short, no M inside non-K/R position + assert qc.synthesis_difficulty in ("LOW", "MODERATE") + + def test_many_flags_high_difficulty(self): + # Long peptide with Met, many Lys/Arg, deamidation sites + qc = check_sequence("s", "GIGKFLKSAKKFGKAFVPEIMNS") + # 23-mer with Met → HIGH difficulty expected + assert qc.synthesis_difficulty == "HIGH" + + def test_no_flags_low(self): + # KKK: charge ~2.8 at pH7.4, ≤2 trypsin sites, no Met/Cys/aggregation/deamidation + qc = check_sequence("s", "KKK") + assert qc.synthesis_difficulty == "LOW" + + +# --------------------------------------------------------------------------- +# check_panel (batch function) +# --------------------------------------------------------------------------- + +class TestCheckPanel: + def test_returns_list_of_qc(self): + panel = [ + {"candidate_id": "A", "sequence": "KRLFKKIGSALKFL"}, + {"candidate_id": "B", "sequence": "RRWQWRMKKLG"}, + ] + results = check_panel(panel) + assert len(results) == 2 + assert all(isinstance(r, SynthQC) for r in results) + + def test_mu_h_map_applied(self): + panel = [{"candidate_id": "A", "sequence": "AAAAAA"}] + mu_h_map = {"A": 0.85} + results = check_panel(panel, mu_h_map=mu_h_map) + assert results[0].mu_h == 0.85 + + def test_empty_panel(self): + assert check_panel([]) == [] + + def test_preserves_candidate_id(self): + panel = [{"candidate_id": "SEED-003_VAR_016", "sequence": "RRWNWRMKKMG"}] + results = check_panel(panel) + assert results[0].candidate_id == "SEED-003_VAR_016" + + def test_missing_mu_h_defaults_to_zero(self): + panel = [{"candidate_id": "A", "sequence": "AAAAAA"}] + results = check_panel(panel) + assert results[0].mu_h == 0.0