diff --git a/CHANGELOG.md b/CHANGELOG.md index cbb4568b..52e0b705 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,8 +9,22 @@ All notable changes to OpenAMP Foundry are documented here. **Status:** Pipeline frozen for synthesis batch ordering. All changes below were quality improvements made before committing the ~$10k wet-lab synthesis budget. -### Latest fixes (PRs #61–69) - +### Latest fixes (PRs #61–70) + +- **PR #70** — Windowed hydrophobic moment (window=11, Eisenberg standard) + anionic charge guard: + `max_windowed_hydrophobic_moment()` added to `features/physchem.py`; `max_hydrophobic_moment` + feature now in `compute_features()` output; `activity_likeness_score()` uses + `max(hydrophobic_moment, max_hydrophobic_moment)` so long-sequence AMPs (magainin-2: 1.51× + improvement, cecropin-A: 2.90×) get credit for their best helical window rather than a + diluted full-sequence average; anionic guard returns 0.0 for peptides with + `charge_density_ph74 < 0.0` (bacterial membrane repulsion; corrects EEIEIEIEIEIEIEE + receiving 0.42 erroneously); fixed Python eager-evaluation bug in feature dict fallback + (`dict.get(key, features["k"])` → `if key in features` conditional); 23 new tests covering + windowed mu_h behaviour, max_hydrophobic_moment in compute_features(), anionic guard cases; + AUROC pipeline.yaml 0.8047→0.8348, phase3.yaml 0.7846→0.8126; regression sentinel updated + (≥0.79); config-identity sentinel updated (<0.82, sits between phase3=0.8126 and + pipeline=0.8348); METHODS.md and DISCOVERY_PREDICTION.md AUROC values updated; + 1223 tests passing - **PR #69** — synthesis_readiness_report.md updated for current 20-candidate pilot panel (AUROC 0.8047; 43 AMPs; 6 seed families; all output files regenerated: presynth_qc_report.md, diversity_report.md, gold_standard_calibration.md); external predictor gate marked PENDING for diff --git a/docs/DISCOVERY_PREDICTION.md b/docs/DISCOVERY_PREDICTION.md index 1458448e..e7bed62e 100644 --- a/docs/DISCOVERY_PREDICTION.md +++ b/docs/DISCOVERY_PREDICTION.md @@ -3,7 +3,7 @@ **Pipeline:** OpenAMP-Foundry v0.1.0 **Date:** 2026-06-28 (updated 2026-06-28) **Status:** Pre-synthesis scientific assessment — for expert review before ordering -**Completed improvements:** Serum stability scoring (PR #31/#32), Family diversity cap (PR #31), Reference set expansion 44→73 sequences (PR #33), Net charge pH 7.4 (PR #34), Helix propensity (PR #35), C-amidation flag (PR #36), Novelty bonus in pilot priority (PR #37), SEED-006 Mastoparan-X (PR #38), Charge×amphipathicity cross-term (PR #39), Amphipathicity weight + helix_bonus (PR #40), SEED-007 Bombolitin II + SEED-008 Puroindoline-a (PR #41), N-terminal acetylation flag + D-amino Wave 2 guidance (PR #42), Full synthesis risk QC — QG/QS deamidation, DG/DS isomerization, Trp photolability (PR #44), Vendor-ready synthesis order generator (PR #45), Diversity-aware pilot panel selection similarity-threshold=0.75 (PR #46), Selectivity proxy + HIGH_CYTOTOX_RISK flag — charge/GRAVY-based mammalian cytotoxicity risk detector (PR #47), selectivity_proxy routed into pilot_priority formula — low-cytotox-risk candidates gently demoted (PR #48), Elastase resistance (HNE 3-protease stability model) + aggregation propensity scoring (synthesis feasibility penalty) (PR #49), Aggregation-safe mutation generation + balanced K/R charge variants + SynthQC continuous aggregation score (PR #50), Proline synthesis penalty + helix bonus enhancement — AUROC 0.8138 → 0.8164 (PR #51), Safety pH74 charge consistency + stronger cytotoxicity penalty — SEED-004 net selectivity swing increased 0.035→0.055 (PR #53), AUPRC metric + SEED-009 (Bac2A proline-rich) + SEED-010 (KR-12 LL-37 fragment) (PR #58), ATCC strain corrections + preregistration assay protocol + negative result schema (PR #60), **CRITICAL FIX: SEED-008 Trp-rich scaffold reinstated in synthesis pool (was excluded by Boman W scale artifact); all 10 seeds generated; SEED-007/009 newly included; SEED-004 correctly excluded by safety gate; 6 mechanism-diverse scaffolds confirmed in synthesis pool** (PR #61), Phase3.yaml AUROC documented (0.7936→0.7890 after PR #65) (PR #63), Stale threshold/reference-count docs corrected (PR #64), **Trp-weighted aromatic bonus (1.5× Trp vs Phe/Tyr; Wimley-White interfacial mechanism); safety abs() bug fix (negative charge does not cause hemolysis); Eisenberg scale comment fix; pool SEED-006: 11→10, SEED-009: 18→19; AUROC 0.8164→0.8086** (PR #65), **Duplicate validation entry removed (REF-GIG-001 = magainin-2 = REF-MAG-001 counted twice); corrected AUROC 0.8086→0.8047 (pipeline), 0.7890→0.7846 (phase3); 43-AMP + 44-background set (n=87)** (PR #66) +**Completed improvements:** Serum stability scoring (PR #31/#32), Family diversity cap (PR #31), Reference set expansion 44→73 sequences (PR #33), Net charge pH 7.4 (PR #34), Helix propensity (PR #35), C-amidation flag (PR #36), Novelty bonus in pilot priority (PR #37), SEED-006 Mastoparan-X (PR #38), Charge×amphipathicity cross-term (PR #39), Amphipathicity weight + helix_bonus (PR #40), SEED-007 Bombolitin II + SEED-008 Puroindoline-a (PR #41), N-terminal acetylation flag + D-amino Wave 2 guidance (PR #42), Full synthesis risk QC — QG/QS deamidation, DG/DS isomerization, Trp photolability (PR #44), Vendor-ready synthesis order generator (PR #45), Diversity-aware pilot panel selection similarity-threshold=0.75 (PR #46), Selectivity proxy + HIGH_CYTOTOX_RISK flag — charge/GRAVY-based mammalian cytotoxicity risk detector (PR #47), selectivity_proxy routed into pilot_priority formula — low-cytotox-risk candidates gently demoted (PR #48), Elastase resistance (HNE 3-protease stability model) + aggregation propensity scoring (synthesis feasibility penalty) (PR #49), Aggregation-safe mutation generation + balanced K/R charge variants + SynthQC continuous aggregation score (PR #50), Proline synthesis penalty + helix bonus enhancement — AUROC 0.8138 → 0.8164 (PR #51), Safety pH74 charge consistency + stronger cytotoxicity penalty — SEED-004 net selectivity swing increased 0.035→0.055 (PR #53), AUPRC metric + SEED-009 (Bac2A proline-rich) + SEED-010 (KR-12 LL-37 fragment) (PR #58), ATCC strain corrections + preregistration assay protocol + negative result schema (PR #60), **CRITICAL FIX: SEED-008 Trp-rich scaffold reinstated in synthesis pool (was excluded by Boman W scale artifact); all 10 seeds generated; SEED-007/009 newly included; SEED-004 correctly excluded by safety gate; 6 mechanism-diverse scaffolds confirmed in synthesis pool** (PR #61), Phase3.yaml AUROC documented (0.7936→0.7890 after PR #65) (PR #63), Stale threshold/reference-count docs corrected (PR #64), **Trp-weighted aromatic bonus (1.5× Trp vs Phe/Tyr; Wimley-White interfacial mechanism); safety abs() bug fix (negative charge does not cause hemolysis); Eisenberg scale comment fix; pool SEED-006: 11→10, SEED-009: 18→19; AUROC 0.8164→0.8086** (PR #65), **Duplicate validation entry removed (REF-GIG-001 = magainin-2 = REF-MAG-001 counted twice); corrected AUROC 0.8086→0.8047 (pipeline), 0.7890→0.7846 (phase3); 43-AMP + 44-background set (n=87)** (PR #66), **Windowed hydrophobic moment (window=11, Eisenberg standard) replaces full-sequence mu_h for activity scoring; anionic guard (charge_density_ph74 < 0.0 → score 0.0); Python eager-eval bug fix in feature dict fallback; 23 new tests; AUROC 0.8047→0.8348 (pipeline), 0.7846→0.8126 (phase3)** (PR #70) --- @@ -243,7 +243,7 @@ Fmoc SPPS with acetonitrile/water gradient; verify MALDI-TOF pre-assay. Updated from ~55–65% due to addition of proven scaffold families. Basis: -- Pipeline AUROC = 0.8047 (95% CI: 0.71–0.89) vs composition-matched UniProt decoys (PR #66 dedup fix; 43 unique AMPs, n=87 total; previous 0.8086 pre-dedup, 0.8164 pre-Trp-bonus) +- Pipeline AUROC = 0.8348 (95% CI: TBD; PR #70 windowed mu_h + anionic guard; 43 unique AMPs, n=87 total; previous 0.8047 PR #66, 0.8086 pre-dedup, 0.8164 pre-Trp-bonus) - Pipeline AUPRC = 0.8443 (+0.2943 above random baseline 0.4943; PR curve emphasises precision at selection operating point) - Recall@20 = 43% on internal benchmark (positives recovered in top 20 ranked candidates) - SEED-003 family (up to 4/20): RRWQWRMKKLG is curated known AMP → variants ~65–75% hit rate @@ -607,7 +607,7 @@ Executing all four actions on the best Wave 1 hits would push the combined proba | Stage | Gate | Original | After PRs #31–38 | After PRs #39–42 | After PRs #43–47 | After PRs #48–53 (current) | Primary limiting factor | |-------|------|----------|-----------------|-----------------|-----------------|--------------------------|------------------------| | 0 | Synthesis success | ~90% | ~90% | ~88% | **~89%** ✓ | **~90%** ✓ (agg model + agg-safe gen + pro penalty + pH74 charge) | SEED-008 W-rich; all aggregation/synthesis risks modelled | -| 1 | MIC ≤ 32 μg/mL | ~55–65% | ~55–65% | ~60–70% | **~60–70%** | **~61–71%** ✓ (AUROC 0.8047) | AUROC 0.8047 (AUPRC 0.8443); 6 scaffold families confirmed | +| 1 | MIC ≤ 32 μg/mL | ~55–65% | ~55–65% | ~60–70% | **~60–70%** | **~62–72%** ✓ (AUROC 0.8348) | AUROC 0.8348 (AUPRC 0.8443); 6 scaffold families confirmed | | 2 | TI > 10 (selectivity) | ~35–50% | ~35–50% | ~38–52% | **~40–55%** ✓ | **~42–57%** ✓ (stronger SEED-004 demotion) | sel_proxy doubled penalty for HIGH_CYTOTOX_RISK tier | | 3 | t½ > 2 h (serum) | ~10–20% | ~25–40% | ~28–42% | **~28–42%** | **~30–46%** ✓ (3-protease model + short/Trp-rich model correction; pilot-panel data) | SEED-003/008 may outperform model score; all seeds borderline — early serum assay recommended | | 4 | Scaffold novelty | ~10–15% | ~18–25% | ~25–35% | **~26–36%** ✓ | **~26–36%** | Diversity filter removes cross-seed near-dups | @@ -683,7 +683,7 @@ The corrected estimate is more conservative because: ### 2. Benchmark Limitations -The AUROC 0.8047 (phase3 synthesis gate: 0.7846) is measured on a small 43+44 demo dataset, not validated against the full +The AUROC 0.8348 (phase3 synthesis gate: 0.8126) is measured on a small 43+44 demo dataset, not validated against the full APD3 (> 3,000 AMPs), DRAMP v3.0 (> 19,000 entries), or ESCAPE benchmark (> 80,000 peptides from 27 repositories). This may overestimate discriminative power. @@ -739,7 +739,7 @@ than the 10 originally planned. The path to 50%+ requires wet-lab data integrati ## Confidence Calibration This assessment is based on: -- Internal benchmark AUROC = 0.8047 (n=87, bootstrap n=2000; phase3 synthesis gate = 0.7846; PR #66 dedup fix) +- Internal benchmark AUROC = 0.8348 (n=87, bootstrap n=2000; phase3 synthesis gate = 0.8126; PR #70 windowed mu_h + anionic guard) - Literature hit rates for physchem AMP prediction (Loose et al. 2006; Tossi et al. 2002) - Published serum stability data for short cationic peptides (Hilpert et al. 2006) - D-amino acid t½ extension data (Wade et al. 1990, PNAS) diff --git a/docs/METHODS.md b/docs/METHODS.md index 7feccf3b..85a52ad1 100644 --- a/docs/METHODS.md +++ b/docs/METHODS.md @@ -102,20 +102,31 @@ seven terms; weights encode biochemical priors only (not data-optimised). Score **Formula (Python source: `scoring/activity.py`):** ``` +# Anionic guard (PR #70): if charge_density_ph74 < 0.0, return 0.0 immediately. +# Anionic peptides are electrostatically repelled by bacterial membranes. + +μH_eff = max(hydrophobic_moment, max_hydrophobic_moment) # windowed mu_h (PR #70) + score = 0.24 × length_score + 0.27 × clamp01((charge_density_ph74 + 0.05) / 0.55) + 0.17 × (1 − min(|hydrophobic − 0.45| / 0.45, 1.0)) - + min(aromatic / 0.20, 1.0) × 0.10 # aromatic_bonus - + clamp01(μH / 0.8) × 0.14 # amphipathicity_score - + clamp01((Pα − 1.0) / 0.20) × 0.03 # helix_bonus - + clamp01(charge_density_ph74 × μH / 0.15) × 0.02 # cross_bonus + + min(aromatic / 0.20, 1.0) × 0.10 # aromatic_bonus + + clamp01(μH_eff / 0.8) × 0.14 # amphipathicity_score (windowed) + + clamp01((Pα − 1.0) / 0.20) × 0.03 # helix_bonus + + clamp01(charge_density_ph74 × μH_eff / 0.15) × 0.02 # cross_bonus # Maximum achievable: 0.24+0.27+0.17+0.10+0.14+0.03+0.02 = 0.97 (arithmetic ceiling; # no explicit min() cap — the code returns clamp01(score) which clips to 1.0) ``` +`μH_eff` = `max(hydrophobic_moment, max_hydrophobic_moment)`. `max_hydrophobic_moment` is the +best μH over any 11-residue window (Eisenberg standard). For sequences ≤ 11 AA this equals +the full-sequence μH. For longer sequences it captures the most amphipathic helical segment, +which better reflects the membrane-active region. (Added PR #70.) + `charge_density_ph74` is the net peptide charge at pH 7.4 (Henderson-Hasselbalch) divided by sequence length; the code falls back to `charge_density` (integer net-charge proxy / length) -for callers that pre-date the pH-7.4 charge model. +for callers that pre-date the pH-7.4 charge model. Lookup uses explicit `if key in features` +conditional (not `dict.get()`) to avoid Python's eager evaluation of the fallback argument. ### 4.1b Boman Activity Score (Second Independent Scorer) @@ -316,29 +327,26 @@ Phase 2 benchmarks verified that the pipeline: - Produces stable rankings under repeated runs (reproducibility) - Shows performance degradation when key scoring dimensions are ablated -**Retrospective AUROC — pipeline.yaml (v0.7.x):** `0.8047` (95% CI: 0.7072–0.8922, n=2000 -bootstrap). Positive-vs-negative separation on the 43-AMP + 44-background benchmark set using -the full ensemble scorer with pipeline.yaml weights. Historical progression: 0.7926 → 0.8138 → +**Retrospective AUROC — pipeline.yaml (v0.8.x):** `0.8348` (95% CI: TBD, n=2000 bootstrap). +Positive-vs-negative separation on the 43-AMP + 44-background benchmark set using the full +ensemble scorer with pipeline.yaml weights. Historical progression: 0.7926 → 0.8138 → 0.8164 (PRs #48–54) → 0.8086 (PR #65: Trp-weighted aromatic bonus, safety abs() fix) → -0.8047 (PR #66: removed duplicate GIGKFLHSAKKFGKAFVGEIMNS / REF-GIG-001 from validation set; -43 unique AMPs used). The minor decrease from 0.8086 reflects that the removed entry (magainin-2) -is a high-ranking AMP; fewer true positives at the top slightly reduce AUROC. Scientifically -more accurate than the pre-dedup estimate. +0.8047 (PR #66: removed duplicate REF-GIG-001; 43 unique AMPs) → +0.8348 (PR #70: windowed mu_h + anionic guard; see below). -**Retrospective AUROC — phase3.yaml (synthesis gate, v0.7.x):** `0.7846` (95% CI: 0.6829–0.8747, +**Retrospective AUROC — phase3.yaml (synthesis gate, v0.8.x):** `0.8126` (95% CI: TBD, n=2000 bootstrap). Phase3 uses re-weighted ensemble scores (activity=0.35, safety=0.30, synthesis=0.20, novelty=0.15 vs pipeline.yaml activity=0.40, safety=0.25, synthesis=0.15, -novelty=0.20) and a stricter safety gate (max_safety_risk=0.40). The AUROC is lower than -pipeline.yaml because the higher safety weight down-ranks some literature AMPs that have -hemolysis risk, which is scientifically correct behaviour for a synthesis gate. Interpretation: -**STRONG** (AUROC > 0.70). AUPRC = 0.8240 (+0.3297 above random baseline of 0.4943). -Recall@10 = 0.23, Recall@20 = 0.44, Recall@43 = 0.70. - -> **Important:** The synthesis selection gate uses phase3.yaml. The phase3 AUROC (0.7846) is the -> operationally relevant benchmark. pipeline.yaml AUROC (0.8047) is the full-ensemble reference. -> Both point estimates exceed the AUROC > 0.70 synthesis gate. Note that with n=87 sequences the -> 95% CI spans the gate boundary; synthesis decisions are made on point estimates, which is standard -> practice at this sample size — the CI reflects sampling uncertainty, not model unreliability. +novelty=0.20) and a stricter safety gate (max_safety_risk=0.40). Historical: 0.7846 (PR #66) +→ 0.8126 (PR #70). The AUROC is lower than pipeline.yaml because the higher safety weight +down-ranks some literature AMPs that have hemolysis risk, which is scientifically correct +behaviour for a synthesis gate. Interpretation: **STRONG** (AUROC > 0.70). + +> **Important:** The synthesis selection gate uses phase3.yaml. The phase3 AUROC (0.8126) is the +> operationally relevant benchmark. pipeline.yaml AUROC (0.8348) is the full-ensemble reference. +> Both point estimates comfortably exceed the AUROC > 0.70 synthesis gate. Note that with n=87 +> sequences the 95% CI spans ~±0.09; synthesis decisions are made on point estimates, which is +> standard practice at this sample size — the CI reflects sampling uncertainty, not model unreliability. **AUPRC (v0.7.x):** `0.8443` for pipeline.yaml (+0.2943 above the random baseline of 0.4943). Area Under Precision-Recall Curve is computed alongside AUROC and reported by `make validate-scoring`. @@ -352,7 +360,7 @@ pessimistic tie-breaking (negatives sort first on tied scores) to match sklearn `make validate-scoring-phase3` (synthesis gate weights, phase3.yaml) or `make validate-scoring-strict` for the composition-matched (scrambled-decoy) variant. -**Key scoring changes in v0.2.x–v0.7.x (PRs #48–65):** +**Key scoring changes in v0.2.x–v0.8.x (PRs #48–70):** | PR | Change | AUROC impact | |----|--------|-------------| @@ -365,6 +373,22 @@ pessimistic tie-breaking (negatives sort first on tied scores) to match sklearn | #58 | AUPRC added; SEED-009 (Bac2A) + SEED-010 (KR-12) added | +2 scaffolds | | #65 | Trp-weighted aromatic bonus (1.5× vs Phe/Tyr); safety abs() fix | 0.8164→0.8086 | | #66 | Remove duplicate REF-GIG-001 (magainin-2 counted twice in validation set) | 0.8086→0.8047 | +| #70 | Windowed mu_h (window=11, Eisenberg standard); anionic charge guard | 0.8047→0.8348 (pipeline); 0.7846→0.8126 (phase3) | + +**PR #70 detail — windowed hydrophobic moment:** For sequences > 11 residues, the full-sequence +μH is diluted by non-helical termini and linker regions. `max_windowed_hydrophobic_moment()` +computes the maximum μH over all 11-residue windows (the Eisenberg standard helical period). +For magainin-2 (23 AA): full-seq μH = 0.4548 vs windowed μH = 0.6862 (1.51×). For cecropin-A +(37 AA): full-seq μH = 0.2267 vs windowed μH = 0.6567 (2.90×). Short sequences (≤ 11 AA) +are unaffected: windowed == full-seq. + +**PR #70 detail — anionic charge guard:** Anionic peptides (charge_density_ph74 < 0.0) are +electrostatically repelled by bacterial membranes (phosphatidylglycerol, cardiolipin surface) +and cannot operate via the canonical cationic AMP mechanism. The guard returns 0.0 immediately, +preventing partial credit from hydrophobic/length terms. Boundary: exactly 0.0 charge density +is not blocked (reflects neutral peptides rather than anionic). Fixed a Python eagerness bug +in feature dict lookup (`features.get(key, features["fallback"])` evaluates fallback eagerly; +replaced with explicit `if "key" in features` conditional). **Limitation:** Benchmarks use a small, curated demo dataset (43 AMPs + 44 background). They do not validate against large independent AMP databases. Real retrospective validation against diff --git a/src/openamp_foundry/features/physchem.py b/src/openamp_foundry/features/physchem.py index 70a78a46..826d69cc 100644 --- a/src/openamp_foundry/features/physchem.py +++ b/src/openamp_foundry/features/physchem.py @@ -139,6 +139,51 @@ def hydrophobic_moment(sequence: str, angle_deg: float = 100.0) -> float: return round(moment, 4) +def max_windowed_hydrophobic_moment( + sequence: str, + window: int = 11, + angle_deg: float = 100.0, +) -> float: + """Maximum hydrophobic moment over any contiguous window of `window` residues. + + Returns a value in [0, ∞); NOT bounded at 1.0 — individual residue Eisenberg + values can be large for tightly packed hydrophobic windows. + + For short sequences (len <= window), this equals the full-sequence mu_h because + only one window exists. For longer sequences, this captures the most concentrated + amphipathic helical segment (Eisenberg 1984, window=11) rather than averaging over + the entire chain, which dilutes the signal near non-helical termini or linkers. + + NOTE: `max_windowed_hydrophobic_moment` is NOT guaranteed to be >= `hydrophobic_moment` + (the full-sequence value). For non-amphipathic or uniform sequences longer than the + window, the fixed-window value (normalised by `window`) can be LESS than the + full-sequence average (normalised by `len(sequence)`). `activity_likeness_score()` + handles this correctly via `max(hydrophobic_moment, max_hydrophobic_moment)`; + direct callers that read both keys from the feature dict must not assume otherwise. + + This value is stored as `max_hydrophobic_moment` in the feature dict and used by + `activity_likeness_score()` to improve discrimination for longer helical AMPs such + as magainin-2, cecropin-A, and LL-37 class peptides. + """ + if not sequence: + return 0.0 + w = min(window, len(sequence)) + angle_rad = math.radians(angle_deg) + best = 0.0 + for start in range(len(sequence) - w + 1): + sub = sequence[start : start + w] + sin_sum = cos_sum = 0.0 + for i, aa in enumerate(sub): + 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) / w + if moment > best: + best = moment + return round(best, 4) + + def interior_protease_sites(sequence: str, cleavage_set: set[str]) -> int: """Count interior residues susceptible to a protease (excludes C-terminal position). @@ -251,6 +296,7 @@ def compute_features(sequence: str) -> dict[str, float | int | dict[str, int]]: pro_fraction = counts.get("P", 0) / length if length else 0.0 repeat_run = longest_repeat_run(sequence) mu_h = hydrophobic_moment(sequence) + max_mu_h = max_windowed_hydrophobic_moment(sequence) n_trypsin = interior_protease_sites(sequence, TRYPSIN_SITES) n_chymotrypsin = interior_protease_sites(sequence, CHYMOTRYPSIN_SITES) n_elastase = interior_protease_sites(sequence, ELASTASE_SITES) @@ -276,6 +322,7 @@ def compute_features(sequence: str) -> dict[str, float | int | dict[str, int]]: "proline_fraction": round(pro_fraction, 4), "longest_repeat_run": repeat_run, "hydrophobic_moment": mu_h, + "max_hydrophobic_moment": max_mu_h, "helix_propensity": helix_pa, "boman_index": boman_index(sequence), "gravy": gravy, diff --git a/src/openamp_foundry/scoring/activity.py b/src/openamp_foundry/scoring/activity.py index d4fbbaf0..6e05d23f 100644 --- a/src/openamp_foundry/scoring/activity.py +++ b/src/openamp_foundry/scoring/activity.py @@ -23,10 +23,32 @@ def activity_likeness_score(features: dict) -> float: length = features["length"] # Prefer physiologically accurate pH-7.4 charge density; fall back to proxy for # manually constructed feature dicts that predate this field. - charge_density = features.get("charge_density_ph74", features["charge_density"]) + # NOTE: dict.get(key, features["fallback"]) evaluates the fallback eagerly in Python; + # use explicit conditional to avoid KeyError when only charge_density_ph74 is present. + if "charge_density_ph74" in features: + charge_density = features["charge_density_ph74"] + else: + charge_density = features["charge_density"] + + # Anionic peptides (net charge < 0 at pH 7.4) are electrostatically repelled by + # bacterial membranes (phosphatidylglycerol, cardiolipin surface) and cannot operate + # via the canonical cationic AMP mechanism modelled by this scorer. + # Return 0.0 immediately rather than awarding partial credit for hydrophobic/length terms. + if charge_density < 0.0: + return 0.0 + hydrophobic = features["hydrophobic_fraction"] aromatic = features["aromatic_fraction"] - mu_h = features.get("hydrophobic_moment", 0.0) + # Use the maximum windowed mu_h (window=11, Eisenberg standard) when available. + # For short sequences (≤11 AA) this equals the full-sequence mu_h. + # For longer sequences it captures the best amphipathic helical segment rather than + # diluting the signal by averaging over non-helical termini or linkers. + # The full-sequence hydrophobic_moment is still used by safety_score() (calibrated + # for its own threshold) — only the activity amphipathicity term changes here. + mu_h = max( + features.get("hydrophobic_moment", 0.0), + features.get("max_hydrophobic_moment", 0.0), + ) # Chou-Fasman mean Pα: range ~0.57 (all Pro/Gly) to ~1.51 (all Glu) # AMP-relevant range: 1.0 (indifferent) to ~1.20 (good helical sequence) helix_pa = features.get("helix_propensity", 1.0) diff --git a/tests/test_physchem_amphipathicity.py b/tests/test_physchem_amphipathicity.py index 7439bba7..4b7be925 100644 --- a/tests/test_physchem_amphipathicity.py +++ b/tests/test_physchem_amphipathicity.py @@ -1,8 +1,13 @@ """Tests for hydrophobic moment (amphipathicity) feature.""" from __future__ import annotations +import pytest -from openamp_foundry.features.physchem import compute_features, hydrophobic_moment +from openamp_foundry.features.physchem import ( + compute_features, + hydrophobic_moment, + max_windowed_hydrophobic_moment, +) class TestHydrophobicMoment: @@ -62,3 +67,118 @@ def test_all_canonical_amino_acids(self): features = compute_features(seq) assert "hydrophobic_moment" in features assert features["hydrophobic_moment"] >= 0.0 + + +class TestMaxWindowedHydrophobicMoment: + """Tests for the Eisenberg-standard windowed mu_h (window=11).""" + + def test_empty_sequence_returns_zero(self): + assert max_windowed_hydrophobic_moment("") == 0.0 + + def test_short_sequence_equals_full_seq_mu_h(self): + # For sequences <= 11 AA, windowed == full-sequence (only one window) + seq = "KWKLFKK" # 7 residues < 11 + windowed = max_windowed_hydrophobic_moment(seq, window=11) + full = hydrophobic_moment(seq) + assert windowed == pytest.approx(full, abs=1e-3) + + def test_long_sequence_windowed_ge_full_seq(self): + # For sequences > 11 AA, windowed >= full-sequence by definition + # (best 11-residue window >= average over all residues) + seq = "GIGKFLHSAKKFGKAFVGEIMNS" # magainin-2, 23 residues + windowed = max_windowed_hydrophobic_moment(seq, window=11) + full = hydrophobic_moment(seq) + assert windowed >= full - 1e-4, ( + f"Windowed mu_h ({windowed:.4f}) should be >= full-seq mu_h ({full:.4f})" + ) + + def test_magainin2_windowed_substantially_higher_than_full(self): + # Magainin-2: windowed mu_h should be meaningfully higher than full-seq + # (helical segment is concentrated; full-seq is diluted by terminal coil) + magainin2 = "GIGKFLHSAKKFGKAFVGEIMNS" + windowed = max_windowed_hydrophobic_moment(magainin2, window=11) + full = hydrophobic_moment(magainin2) + # Empirically: windowed ~0.68, full ~0.45 → ratio ~1.5× + assert windowed > full * 1.2, ( + f"Magainin-2 windowed ({windowed:.4f}) should be >1.2× full-seq ({full:.4f})" + ) + + def test_returns_float_rounded_to_4dp(self): + result = max_windowed_hydrophobic_moment("KWKLFKKIGAVLKVL") + assert isinstance(result, float) + assert result == round(result, 4) + + def test_single_residue_equals_full_moment(self): + result = max_windowed_hydrophobic_moment("K", window=11) + full = hydrophobic_moment("K") + assert result == pytest.approx(full, abs=1e-3) + + def test_window_equal_to_sequence_length(self): + seq = "KWKLFKKIGAV" # exactly 11 residues + windowed = max_windowed_hydrophobic_moment(seq, window=11) + full = hydrophobic_moment(seq) + assert windowed == pytest.approx(full, abs=1e-3) + + def test_nonnegative(self): + for seq in ["", "A", "KWKLFKK", "GIGKFLHSAKKFGKAFVGEIMNS", "EEIEIEIEIEIEIEE"]: + result = max_windowed_hydrophobic_moment(seq) + assert result >= 0.0, f"max_windowed_hydrophobic_moment({seq!r}) returned {result}" + + def test_custom_window(self): + seq = "KWKLFKKIGAVLKVL" + w7 = max_windowed_hydrophobic_moment(seq, window=7) + w11 = max_windowed_hydrophobic_moment(seq, window=11) + # Both must be non-negative; neither is bounded at 1.0 (Eisenberg values can + # be large for tightly packed hydrophobic windows — clamping happens in the scorer) + assert w7 >= 0.0 + assert w11 >= 0.0 + + +class TestComputeFeaturesMaxMuH: + """Tests that compute_features() includes max_hydrophobic_moment.""" + + def test_max_hydrophobic_moment_key_present(self): + features = compute_features("KWKLFKKIGAVLKVL") + assert "max_hydrophobic_moment" in features, ( + "compute_features() must return 'max_hydrophobic_moment'" + ) + + def test_max_hydrophobic_moment_is_float(self): + features = compute_features("KWKLFKKIGAVLKVL") + assert isinstance(features["max_hydrophobic_moment"], float) + + def test_max_hydrophobic_moment_nonnegative(self): + for seq in ["", "A", "KWKLFKK", "GIGKFLHSAKKFGKAFVGEIMNS"]: + features = compute_features(seq) + assert features["max_hydrophobic_moment"] >= 0.0 + + def test_max_hydrophobic_moment_ge_hydrophobic_moment_for_known_amps(self): + # For known helical AMPs with a concentrated amphipathic segment, the windowed + # value should be >= the full-seq value (good amphipathic signal). + # NOTE: this invariant does NOT hold for non-amphipathic or uniform sequences + # (see test_windowed_can_be_less_than_full_for_uniform_sequences below). + for seq in ["KWKLFKKIGAVLKVL", "GIGKFLHSAKKFGKAFVGEIMNS", "LLGDFFRKSKEKIGKEFKRIVQRIKDFLRNLVPRTES"]: + features = compute_features(seq) + assert features["max_hydrophobic_moment"] >= features["hydrophobic_moment"] - 1e-4, ( + f"max_mu_h ({features['max_hydrophobic_moment']:.4f}) < mu_h " + f"({features['hydrophobic_moment']:.4f}) for {seq}" + ) + + def test_windowed_can_be_less_than_full_for_uniform_sequences(self): + # For non-amphipathic sequences (e.g. all-Ala, all-Glu), the windowed mu_h + # (normalised by window=11) can be LESS than full-seq mu_h (normalised by len>11). + # activity_likeness_score() handles this correctly via max(); direct callers + # must not assume max_hydrophobic_moment >= hydrophobic_moment. + seq = "AAAAAAAAAAAA" # 12-A, uniform → no concentrated amphipathic segment + windowed = max_windowed_hydrophobic_moment(seq, window=11) + full = hydrophobic_moment(seq) + assert windowed < full, ( + f"For uniform sequence {seq!r}, expected windowed ({windowed:.4f}) < " + f"full-seq ({full:.4f}) to document the non-invariant" + ) + + def test_short_seq_max_equals_full(self): + features = compute_features("KWKLFKK") # 7 residues, window=11 → only one window + assert features["max_hydrophobic_moment"] == pytest.approx( + features["hydrophobic_moment"], abs=1e-3 + ) diff --git a/tests/test_retrospective.py b/tests/test_retrospective.py index 0480e118..79c40f29 100644 --- a/tests/test_retrospective.py +++ b/tests/test_retrospective.py @@ -144,11 +144,11 @@ def test_phase3_benchmark_passes_gate(self): phase3.yaml shifts more weight to safety (0.30 vs 0.25) and synthesis (0.20 vs 0.15), which down-ranks hemolytic AMPs — correct behaviour for a synthesis gate but lowers raw AUROC slightly vs pipeline.yaml. Gate is still AUROC > 0.70 (STRONG). - Measured: AUROC=0.7846 (95% CI 0.6829-0.8747, n=2000 bootstrap; after PR #66 dedup fix). - Previous values: 0.7936 (before PR #65 Trp bonus), 0.7890 (pre-dedup, post-PR #65). + Measured: AUROC=0.8126 (post PR #70 windowed mu_h + anionic guard). + Previous values: 0.7846 (pre-PR #70), 0.7936 (before PR #65 Trp bonus). - The upper-bound check (< 0.83) is a config-identity sentinel: pipeline.yaml scores - 0.8047 which would breach it, so a silent config fallback is detectable. + The upper-bound check (< 0.82) is a config-identity sentinel: pipeline.yaml scores + 0.8348 which would breach it, so a silent config fallback is detectable. """ from pathlib import Path amp_csv = Path("examples/validation/known_amps.csv") @@ -166,17 +166,17 @@ def test_phase3_benchmark_passes_gate(self): f"phase3.yaml AUROC={result['auroc']:.4f}: synthesis gate does not meet the 0.70 " "threshold. Do not proceed to wet-lab synthesis." ) - # Regression sentinel: measured AUROC=0.7846 (post-dedup); alert if it drops below 0.75 - assert result["auroc"] >= 0.75, ( - f"phase3.yaml AUROC={result['auroc']:.4f} < 0.75: a regression has occurred " - "(baseline 0.7846). Synthesis gate still valid at >0.70 but scoring may be degraded." + # Regression sentinel: measured AUROC=0.8126 (post PR #70 windowed mu_h + anionic guard) + assert result["auroc"] >= 0.79, ( + f"phase3.yaml AUROC={result['auroc']:.4f} < 0.79: a regression has occurred " + "(baseline 0.8126). Synthesis gate still valid at >0.70 but scoring may be degraded." ) - # Config-identity sentinel: pipeline.yaml AUROC=0.8047 >= 0.795 would fail this, - # while phase3.yaml AUROC=0.7846 < 0.795 passes. Threshold sits between the two measured + # Config-identity sentinel: pipeline.yaml AUROC=0.8348 >= 0.82 would fail this, + # while phase3.yaml AUROC=0.8126 < 0.82 passes. Threshold sits between the two measured # point estimates so a silent config fallback is caught. - assert result["auroc"] < 0.795, ( - f"phase3.yaml AUROC={result['auroc']:.4f} >= 0.795: this is higher than expected " - "for phase3.yaml (measured 0.7846). Check that phase3.yaml (not pipeline.yaml) was loaded." + assert result["auroc"] < 0.82, ( + f"phase3.yaml AUROC={result['auroc']:.4f} >= 0.82: this is higher than expected " + "for phase3.yaml (measured 0.8126). Check that phase3.yaml (not pipeline.yaml) was loaded." ) print(f"\n[phase3 benchmark] AUROC={result['auroc']:.4f}: {result['interpretation']}") print(f"AUPRC={result.get('auprc', 'N/A')}, Recall@20={result.get('recall_at_20', 'N/A')}") diff --git a/tests/test_scoring_activity.py b/tests/test_scoring_activity.py index 62cc197b..ef552ab5 100644 --- a/tests/test_scoring_activity.py +++ b/tests/test_scoring_activity.py @@ -124,6 +124,68 @@ def test_length_30_midrange_verifies_normalization_constant(self): assert abs(score - 0.4666) < 0.001 +class TestAnionicGuard: + """Anionic peptides must receive score 0.0 — they are electrostatically repelled + by bacterial membranes and cannot operate via the cationic AMP mechanism.""" + + def test_anionic_peptide_returns_zero(self): + # charge_density_ph74 < 0.0 → guard fires, return 0.0 regardless of other terms + feat = { + "length": 16, "charge_density_ph74": -0.5, + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.10, + "hydrophobic_moment": 0.6, + } + assert activity_likeness_score(feat) == 0.0 + + def test_anionic_guard_uses_ph74_field_preferentially(self): + # charge_density_ph74 field takes precedence over legacy charge_density + feat = { + "length": 18, "charge_density": 0.30, # positive legacy field + "charge_density_ph74": -0.3, # but negative at pH 7.4 + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.0, + "hydrophobic_moment": 0.0, + } + assert activity_likeness_score(feat) == 0.0 + + def test_anionic_guard_polyglutamate(self): + # EEIEIEIEIEIEIEE-like peptide: highly structured but anionic — must be 0 + feat = { + "length": 15, "charge_density_ph74": -0.6, + "hydrophobic_fraction": 0.40, "aromatic_fraction": 0.0, + "hydrophobic_moment": 0.42, "max_hydrophobic_moment": 0.55, + } + assert activity_likeness_score(feat) == 0.0 + + def test_zero_charge_density_not_blocked(self): + # Exactly 0.0 charge density should NOT trigger the guard (guard is < 0.0) + feat = { + "length": 18, "charge_density_ph74": 0.0, + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.0, + "hydrophobic_moment": 0.0, + } + score = activity_likeness_score(feat) + # Should still get some score from length + hydrophobic terms + assert score > 0.0 + + def test_small_negative_charge_blocked(self): + # Even mildly negative (e.g. -0.05) must return 0.0 + feat = { + "length": 18, "charge_density_ph74": -0.05, + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.0, + "hydrophobic_moment": 0.0, + } + assert activity_likeness_score(feat) == 0.0 + + def test_positive_charge_not_blocked(self): + # Positive charge density must pass through normally + feat = { + "length": 18, "charge_density": 0.25, + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.0, + "hydrophobic_moment": 0.0, + } + assert activity_likeness_score(feat) > 0.0 + + class TestChargeComponent: def test_negative_charge_penalised(self): neg = activity_likeness_score(_feat(charge_density=-0.5)) @@ -290,6 +352,46 @@ def test_trp_bonus_active_only_in_sub_saturation_range(self): ) +class TestMaxWindowedMuHInActivity: + """activity_likeness_score uses max(hydrophobic_moment, max_hydrophobic_moment) + so that long-sequence AMPs get credit for their best helical window.""" + + def test_max_hydrophobic_moment_used_when_higher(self): + # max_hydrophobic_moment > hydrophobic_moment → score should reflect the max + feat_with_max = { + "length": 23, "charge_density": 0.25, + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.0, + "hydrophobic_moment": 0.45, + "max_hydrophobic_moment": 0.69, # windowed value higher than full-seq + } + feat_full_only = { + "length": 23, "charge_density": 0.25, + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.0, + "hydrophobic_moment": 0.45, + # no max_hydrophobic_moment key + } + assert activity_likeness_score(feat_with_max) > activity_likeness_score(feat_full_only) + + def test_full_seq_used_when_no_max_key(self): + # Backward compat: if max_hydrophobic_moment absent, full-seq mu_h is used + feat = { + "length": 18, "charge_density": 0.30, + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.0, + "hydrophobic_moment": 0.6, + } + score = activity_likeness_score(feat) + assert score > 0.0 + + def test_missing_both_mu_h_keys_defaults_to_zero(self): + feat = { + "length": 18, "charge_density": 0.30, + "hydrophobic_fraction": 0.45, "aromatic_fraction": 0.0, + } + score = activity_likeness_score(feat) + # Should still work; mu_h defaults to 0 → amphipathicity_score = 0 + assert score >= 0.0 + + class TestAmphipathicityBonus: def test_no_amphipathicity_no_bonus(self): s_0 = activity_likeness_score(_feat(mu_h=0.0)) @@ -297,7 +399,7 @@ def test_no_amphipathicity_no_bonus(self): assert s_4 > s_0 def test_mu_h_0_8_full_bonus(self): - # amphipathicity = clamp01(0.8/0.8) * 0.15 = 0.15 + # amphipathicity = clamp01(0.8/0.8) * 0.14 = 0.14 s_08 = activity_likeness_score(_feat(mu_h=0.8)) s_10 = activity_likeness_score(_feat(mu_h=1.0)) # Both give full bonus (capped)