diff --git a/docs/DISCOVERY_PREDICTION.md b/docs/DISCOVERY_PREDICTION.md index 649f3055..5a23a839 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) +**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) --- @@ -243,8 +243,8 @@ 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.8164 (95% CI: 0.72–0.90) vs composition-matched UniProt decoys -- Pipeline AUPRC = 0.8556 (+0.36 above random baseline 0.50; PR curve emphasises precision at selection operating point) +- Pipeline AUROC = 0.8086 (95% CI: 0.71–0.90) vs composition-matched UniProt decoys (PR #65 update; previous 0.8164 after Trp bonus scientific fix; within bootstrap CI) +- Pipeline AUPRC = 0.8495 (+0.35 above random baseline 0.50; 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 - SEED-007 family (up to 4/20): Bombolitin II (*Bombus pennsylvanicus* bumblebee venom) with @@ -543,7 +543,7 @@ in the 20-member pilot panel, P(≥1 from family) ≈ 1 − (1 − P_candidate)^ | 1 | SEED-008 | FPVTWRWWKWWKG vars | 0.385 | **0.667** | ~15–25% | Highest novelty + proven Trp mechanism; 24 selected (PR #61 reinstated) | | 2 | SEED-007 | IKITTMLKKLG vars | 0.636 | 0.615 | ~14–22% | Bumblebee venom; distinct from wasp mastoparan; 27 selected (highest count) | | 3 | SEED-006 | INWKGIAAMAKKLL vars | 0.667 | 0.643 | ~14–22% | Balanced stability+novelty; mastoparan mechanism data | -| 4 | SEED-009 | RRLPRPPYLPRP vars | 0.80+ | 0.467 | ~12–20% | Proline-rich intracellular (DnaK); orthogonal mechanism to all other seeds; 18 selected | +| 4 | SEED-009 | RRLPRPPYLPRP vars | 0.80+ | 0.467 | ~12–20% | Proline-rich intracellular (DnaK); orthogonal mechanism to all other seeds; 19 selected | | 5 | SEED-003 | RRWQWRMKKLG vars | 0.27 | 0.182 | ~8–16% | High ensemble (0.780-0.795) but poor serum stability; 16 selected | | 6 | SEED-005 | KRLFKKIGSALKFL (seed) | 0.52 | 0.467 | ~8–14% | Cecropin-magainin hybrid; best variant VAR_009 (KRFFKKIGSALKFA, FF motif) aids membrane insertion; 4 selected | @@ -584,7 +584,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.8164) | AUROC 0.8164 (AUPRC computed); 10 scaffold families | +| 1 | MIC ≤ 32 μg/mL | ~55–65% | ~55–65% | ~60–70% | **~60–70%** | **~61–71%** ✓ (AUROC 0.8086) | AUROC 0.8086 (AUPRC 0.8495); 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%** | **~29–44%** ✓ (3-protease model) | Wave 2 D-amino plan machine-readable | | 4 | Scaffold novelty | ~10–15% | ~18–25% | ~25–35% | **~26–36%** ✓ | **~26–36%** | Diversity filter removes cross-seed near-dups | @@ -660,7 +660,7 @@ The corrected estimate is more conservative because: ### 2. Benchmark Limitations -The AUROC 0.8164 is measured on a small 44+44 demo dataset, not validated against the full +The AUROC 0.8086 (phase3 synthesis gate: 0.7890) is measured on a small 44+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. @@ -716,7 +716,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.8164 (n=88, bootstrap n=2000; improved from 0.8138 after helix bonus weight PR #51) +- Internal benchmark AUROC = 0.8086 (n=88, bootstrap n=2000; phase3 synthesis gate = 0.7890; PR #65 Trp-weighted bonus update) - 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) @@ -739,8 +739,8 @@ is the fastest path from the current ~10–18% (calibrated) to the ≥20–30% t justify a larger Wave 2 investment. **Synthesis pool status (PR #61, final pre-wetlab state):** -6 scaffold families confirmed: SEED-003 (16), SEED-005 (4), SEED-006 (11), SEED-007 (27), -SEED-008 (24), SEED-009 (18). Total: 100 selected candidates. Mean safety=0.991, mean +6 scaffold families confirmed: SEED-003 (16), SEED-005 (4), SEED-006 (10), SEED-007 (27), +SEED-008 (24), SEED-009 (19). Total: 100 selected candidates. Mean safety=0.991, mean ensemble=0.805. 20-member pilot panel drawn from these 100 for Wave 1 synthesis. --- diff --git a/docs/METHODS.md b/docs/METHODS.md index 7f9f3b41..8cf735f9 100644 --- a/docs/METHODS.md +++ b/docs/METHODS.md @@ -316,29 +316,29 @@ 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.2.x):** `0.8164` (positive-vs-negative separation on -the 44-AMP + 44-background benchmark set using the full ensemble scorer with pipeline.yaml -weights). Benchmarked after PRs #48–54; the scoring improvements in those PRs moved AUROC from -0.7926 → 0.8138 → 0.8164. Bootstrap 95% CI computed via Wilcoxon-Mann-Whitney (n=2000 -iterations; CI reported in `validate-scoring` output). - -**Retrospective AUROC — phase3.yaml (synthesis gate, v0.6.x):** `0.7936` (95% CI: 0.6963–0.8827, +**Retrospective AUROC — pipeline.yaml (v0.7.x):** `0.8086` (95% CI: 0.7149–0.8957, n=2000 +bootstrap). Positive-vs-negative separation on the 44-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). The minor +decrease from 0.8164 reflects that some random-background decoys contain Trp residues and benefit +from the Trp bonus; this is within bootstrap CI width and the bonus is scientifically correct. + +**Retrospective AUROC — phase3.yaml (synthesis gate, v0.7.x):** `0.7890` (95% CI: 0.6932–0.8784, 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.8333 (+0.3333 above random baseline of 0.50). -Recall@10 = 0.23, Recall@20 = 0.43, Recall@44 = 0.68. +**STRONG** (AUROC > 0.70). AUPRC = 0.8301 (+0.3301 above random baseline of 0.50). +Recall@10 = 0.23, Recall@20 = 0.43, Recall@44 = 0.70. -> **Important:** The synthesis selection gate uses phase3.yaml. The phase3 AUROC (0.7936) is the -> operationally relevant benchmark. pipeline.yaml AUROC (0.8164) is the full-ensemble reference. +> **Important:** The synthesis selection gate uses phase3.yaml. The phase3 AUROC (0.7890) is the +> operationally relevant benchmark. pipeline.yaml AUROC (0.8086) is the full-ensemble reference. > Both point estimates exceed the AUROC > 0.70 synthesis gate. Note that with n=88 sequences the -> 95% CI spans the gate boundary (0.6963–0.8827); synthesis decisions are made on point estimates, -> which is standard practice at this sample size — the CI reflects sampling uncertainty, not model -> unreliability. +> 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. -**AUPRC (v0.5.x):** `0.8556` for pipeline.yaml (+0.3556 above the random baseline of 0.50). +**AUPRC (v0.7.x):** `0.8495` for pipeline.yaml (+0.3495 above the random baseline of 0.50). Area Under Precision-Recall Curve is computed alongside AUROC and reported by `make validate-scoring`. AUPRC is preferable to AUROC for class-imbalanced datasets because it emphasises precision at the operating point actually used for candidate selection. Random baseline = dataset @@ -350,7 +350,7 @@ and avoid inflation to 1.0 for constant-score classifiers. `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.5.x (PRs #48–58):** +**Key scoring changes in v0.2.x–v0.7.x (PRs #48–65):** | PR | Change | AUROC impact | |----|--------|-------------| @@ -361,6 +361,7 @@ and avoid inflation to 1.0 for constant-score classifiers. | #51 | Proline synthesis penalty (−0.10 at >15%); helix bonus 0.01→0.03 | 0.8138→0.8164 | | #53 | Safety uses charge_density_ph74; stronger cytotox_penalty in pilot | pilot ranking | | #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 | **Limitation:** Benchmarks use a small, curated demo dataset (44 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 00b7245c..70a78a46 100644 --- a/src/openamp_foundry/features/physchem.py +++ b/src/openamp_foundry/features/physchem.py @@ -19,8 +19,10 @@ # Hydrophobic residues that drive aggregation in synthetic peptides: Val, Ile, Leu, Met, Phe, Trp AGG_HYDROPHOBIC = set("VILMFW") -# 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. +# Eisenberg (1984) consensus hydrophobicity scale. Raw signed values: R=-2.530 (most hydrophilic) +# to I=+1.380 (most hydrophobic). DO NOT normalize to [0,1] — the hydrophobic moment calculation +# requires signed values so that polar/charged residues (negative) cancel apolar residues (positive). +# Source: Eisenberg et al. (1984) J Mol Biol 179:125-142, Table 2. 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, @@ -71,6 +73,12 @@ def net_charge_at_ph74(sequence: str) -> float: Reference: Bjellqvist et al. (1993) Electrophoresis 14:1023; pKa values for His in peptides are shifted vs free amino acid (6.0) to 6.5 to account for context. + + Note: α-amino N-terminus (+≈0.80 at pH 7.4, pKa ≈ 8.0) and α-carboxyl C-terminus + (−≈1.00 at pH 7.4, pKa ≈ 3.1) are excluded — only side-chain charges are computed. + The net terminal offset is ≈ −0.20 per peptide (~7% underestimation for a 10-mer). + This is the standard convention for relative charge comparisons across candidates + synthesized under the same conditions. """ charge = 0.0 for aa in sequence: diff --git a/src/openamp_foundry/scoring/activity.py b/src/openamp_foundry/scoring/activity.py index 88e80bf4..d4fbbaf0 100644 --- a/src/openamp_foundry/scoring/activity.py +++ b/src/openamp_foundry/scoring/activity.py @@ -40,8 +40,17 @@ def activity_likeness_score(features: dict) -> float: # 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 + # Aromatic residues (F, W, Y) aid membrane insertion. + # Trp is weighted 1.5× vs Phe/Tyr: Trp's indole ring preferentially anchors at the + # lipid-water interface (Wimley-White W=-1.85 kcal/mol; Wimley & White 1996 Nat Struct Biol), + # a mechanism distinct from and more membrane-favorable than Phe/Tyr's hydrophobic stacking. + # This differentiates puroindoline-a/indolicidin-class Trp-rich AMPs from Phe-rich peptides. + _counts = features.get("residue_counts", {}) + _length = features["length"] or 1 + trp_fraction = _counts.get("W", 0) / _length + non_trp_aromatic = max(0.0, aromatic - trp_fraction) + weighted_aromatic = trp_fraction * 1.5 + non_trp_aromatic + aromatic_bonus = min(weighted_aromatic / 0.20, 1.0) * 0.10 # 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 diff --git a/src/openamp_foundry/scoring/boman.py b/src/openamp_foundry/scoring/boman.py index 05562020..2401050a 100644 --- a/src/openamp_foundry/scoring/boman.py +++ b/src/openamp_foundry/scoring/boman.py @@ -22,6 +22,8 @@ """ from __future__ import annotations +import math + # Published amino acid interaction potentials from Boman (2003), Table 1, set 2. # Units: kcal/mol (approximate interaction energies) # Positive → hydrophilic/charged, negative → hydrophobic @@ -86,7 +88,6 @@ def boman_activity_score(sequence: str) -> float: bi = boman_index(sequence) # Sigmoid-like normalization centered at 0, scaled so bi=1.0 → ≈ 0.73 # Using tanh mapping: score = 0.5 * (1 + tanh(bi / 2)) - import math score = 0.5 * (1.0 + math.tanh(bi / 2.0)) return round(max(0.0, min(1.0, score)), 4) diff --git a/src/openamp_foundry/scoring/safety.py b/src/openamp_foundry/scoring/safety.py index dd9572e0..42159742 100644 --- a/src/openamp_foundry/scoring/safety.py +++ b/src/openamp_foundry/scoring/safety.py @@ -17,9 +17,10 @@ def safety_score(features: dict) -> float: typical AMP-like candidates (μH ≈ 0.35–0.50). """ length = features["length"] - # Prefer pH 7.4 charge (same key-preference as activity scorer); abs() applied because - # charge magnitude, not sign, drives hemolysis risk. Fall back for pre-PR-47 callers. - charge_density = abs(features.get("charge_density_ph74", features["charge_density"])) + # Prefer pH 7.4 charge density. Positive charge drives hemolysis (electrostatic membrane + # disruption); negatively charged peptides are repelled from both mammalian and bacterial + # membranes and pose no hemolysis risk. Do NOT use abs() here. + charge_density = features.get("charge_density_ph74", features["charge_density"]) hydrophobic = features["hydrophobic_fraction"] cys = features["cysteine_fraction"] repeat_run = features["longest_repeat_run"] diff --git a/tests/test_retrospective.py b/tests/test_retrospective.py index 1076fa63..51f4a3e0 100644 --- a/tests/test_retrospective.py +++ b/tests/test_retrospective.py @@ -144,7 +144,7 @@ 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.7936 (95% CI 0.6963-0.8827, n=2000 bootstrap). + Measured: AUROC=0.7890 (95% CI 0.6932-0.8784, n=2000 bootstrap; after PR #65 Trp bonus). The upper-bound check (< 0.83) is a config-identity sentinel: pipeline.yaml scores 0.8164 which would breach it, so a silent config fallback is detectable. diff --git a/tests/test_scoring_activity.py b/tests/test_scoring_activity.py index b290c264..62cc197b 100644 --- a/tests/test_scoring_activity.py +++ b/tests/test_scoring_activity.py @@ -4,7 +4,9 @@ - length_score (weight 0.24): 1 - min(|length - 18| / 25, 1.0) - charge_score (weight 0.27): clamp01((charge_density + 0.05) / 0.55) - hydro_score (weight 0.17): 1 - min(|hf - 0.45| / 0.45, 1.0) - - aromatic_bonus (max 0.10): min(aromatic / 0.20, 1.0) * 0.10 + - aromatic_bonus (max 0.10): min(weighted_aromatic / 0.20, 1.0) * 0.10 + where weighted_aromatic = trp_fraction * 1.5 + non_trp_aromatic + (Trp weighted 1.5× vs Phe/Tyr; Wimley-White interfacial insertion mechanism) - amphipathicity (max 0.14): clamp01(mu_h / 0.8) * 0.14 - helix_bonus (max 0.03): clamp01((helix_pa - 1.0) / 0.20) * 0.03 - cross_bonus (max 0.02): clamp01(charge_density * mu_h / 0.15) * 0.02 @@ -32,6 +34,25 @@ def _feat( } +def _feat_with_counts( + length: int = 20, + charge_density: float = 0.30, + hydrophobic: float = 0.45, + aromatic: float = 0.0, + mu_h: float = 0.0, + counts: dict | None = None, +) -> dict: + """Feature dict with residue_counts to exercise the Trp-weighted code path.""" + return { + "length": length, + "charge_density": charge_density, + "hydrophobic_fraction": hydrophobic, + "aromatic_fraction": aromatic, + "hydrophobic_moment": mu_h, + "residue_counts": counts or {}, + } + + class TestClamp01: def test_clamps_below_zero(self): assert clamp01(-5.0) == 0.0 @@ -184,6 +205,91 @@ def test_aromatic_half_gives_half_bonus(self): assert abs(diff - 0.05) < 0.001 +class TestTrpWeightedAromaticBonus: + """Tests for the Trp-weighted aromatic bonus (residue_counts code path). + + Trp is weighted 1.5× vs Phe/Tyr because Trp anchors at the lipid-water + interface (Wimley-White interfacial scale W=−1.85 kcal/mol). + """ + + def test_trp_rich_scores_higher_than_phe_rich_same_aromatic_fraction(self): + """Trp-rich peptide should outscore Phe-rich at identical aromatic fraction.""" + # 2 W out of 20 = 10% aromatic; 2 F out of 20 = 10% aromatic + feat_trp = _feat_with_counts(length=20, aromatic=0.10, counts={"W": 2}) + feat_phe = _feat_with_counts(length=20, aromatic=0.10, counts={"F": 2}) + assert activity_likeness_score(feat_trp) > activity_likeness_score(feat_phe), ( + "Trp-rich peptide should score higher than Phe-rich at the same aromatic fraction" + ) + + def test_trp_bonus_factor_1_5x_quantitative(self): + """With only Trp (no Phe/Tyr): weighted_aromatic = trp_fraction * 1.5. + + trp_fraction=0.10 → weighted=0.15 → bonus = (0.15/0.20)*0.10 = 0.075 + vs. old Phe-based (0.10/0.20)*0.10 = 0.05 → delta = 0.025 + """ + feat_trp = _feat_with_counts(length=20, aromatic=0.10, counts={"W": 2}) + feat_phe = _feat_with_counts(length=20, aromatic=0.10, counts={"F": 2}) + delta = activity_likeness_score(feat_trp) - activity_likeness_score(feat_phe) + # Expected: bonus_trp - bonus_phe = 0.075 - 0.05 = 0.025 + assert abs(delta - 0.025) < 0.002, ( + f"Trp vs Phe aromatic bonus delta should be ~0.025, got {delta:.4f}" + ) + + def test_trp_aromatic_bonus_capped_at_0_10(self): + """Ceiling is 0.10 regardless of Trp fraction — 1.5x just reaches cap sooner.""" + # trp_fraction=0.14 → weighted=0.21 > 0.20 → capped: bonus=1.0*0.10=0.10 + feat_trp_cap = _feat_with_counts(length=14, aromatic=0.14, counts={"W": 2}) + # phe_fraction=0.20 → bonus=1.0*0.10=0.10 (also at cap) + feat_phe_cap = _feat_with_counts(length=20, aromatic=0.20, counts={"F": 4}) + # Both should produce the same aromatic contribution (0.10) + # They may differ on other terms (length), so isolate aromatic by comparing + # to zero-aromatic baseline at same length + feat_trp_zero = _feat_with_counts(length=14, aromatic=0.0, counts={}) + feat_phe_zero = _feat_with_counts(length=20, aromatic=0.0, counts={}) + trp_bonus = activity_likeness_score(feat_trp_cap) - activity_likeness_score(feat_trp_zero) + phe_bonus = activity_likeness_score(feat_phe_cap) - activity_likeness_score(feat_phe_zero) + assert abs(trp_bonus - 0.10) < 0.002, f"Trp cap bonus should be 0.10, got {trp_bonus:.4f}" + assert abs(phe_bonus - 0.10) < 0.002, f"Phe cap bonus should be 0.10, got {phe_bonus:.4f}" + + def test_no_residue_counts_falls_back_to_unweighted(self): + """When residue_counts is absent, Trp fraction defaults to 0 → same as all-Phe.""" + # Without residue_counts: weighted_aromatic = 0*1.5 + aromatic = aromatic (Phe behavior) + feat_without_counts = _feat(length=20, aromatic=0.10) + feat_phe = _feat_with_counts(length=20, aromatic=0.10, counts={"F": 2}) + assert activity_likeness_score(feat_without_counts) == activity_likeness_score(feat_phe) + + def test_mixed_trp_phe_intermediate_score(self): + """Sequence with 1 Trp + 1 Phe (same aromatic fraction) scores between all-W and all-F.""" + feat_trp = _feat_with_counts(length=20, aromatic=0.10, counts={"W": 2}) + feat_phe = _feat_with_counts(length=20, aromatic=0.10, counts={"F": 2}) + feat_mixed = _feat_with_counts(length=20, aromatic=0.10, counts={"W": 1, "F": 1}) + score_trp = activity_likeness_score(feat_trp) + score_phe = activity_likeness_score(feat_phe) + score_mixed = activity_likeness_score(feat_mixed) + assert score_phe < score_mixed < score_trp, ( + f"Mixed Trp/Phe should score between all-Phe ({score_phe:.4f}) " + f"and all-Trp ({score_trp:.4f}), got {score_mixed:.4f}" + ) + + def test_trp_bonus_active_only_in_sub_saturation_range(self): + """The Trp bonus differentiates at MODERATE aromatic fraction (< 13% Trp to reach cap). + + At high aromatic fraction (>= 20%), both Trp and Phe hit the cap and the bonus + is identical. SEED-008 (38% Trp) already saturates; the bonus matters for sequences + with 5-13% Trp where the 1.5x weight pushes them above equivalent Phe fractions. + """ + # Low aromatic fraction: 5% Trp vs 5% Phe — Trp gets higher score + feat_low_trp = _feat_with_counts(length=20, aromatic=0.05, counts={"W": 1}) + feat_low_phe = _feat_with_counts(length=20, aromatic=0.05, counts={"F": 1}) + assert activity_likeness_score(feat_low_trp) > activity_likeness_score(feat_low_phe) + # High aromatic fraction: 25% both → both cap at 0.10, no difference + feat_hi_trp = _feat_with_counts(length=20, aromatic=0.25, counts={"W": 5}) + feat_hi_phe = _feat_with_counts(length=20, aromatic=0.25, counts={"F": 5}) + assert activity_likeness_score(feat_hi_trp) == activity_likeness_score(feat_hi_phe), ( + "High aromatic fraction: both Trp and Phe should saturate the cap at 0.10" + ) + + class TestAmphipathicityBonus: def test_no_amphipathicity_no_bonus(self): s_0 = activity_likeness_score(_feat(mu_h=0.0)) diff --git a/tests/test_scoring_safety.py b/tests/test_scoring_safety.py index e891288a..fbf0001a 100644 --- a/tests/test_scoring_safety.py +++ b/tests/test_scoring_safety.py @@ -3,7 +3,7 @@ safety_score = 1.0 - clamp01(risk), where risk accumulates from: - μH > 0.55: (μH - 0.55) * 1.5 - hydrophobic_fraction > 0.65: (hf - 0.65) * 1.8 - - |charge_density| > 0.55: (cd - 0.55) * 1.2 + - charge_density > 0.55: (cd - 0.55) * 1.2 [positive charge only — abs() removed] - length > 35: +0.25 flat - cysteine_fraction > 0.25: +0.20 flat - longest_repeat_run >= 6: +0.15 flat @@ -114,14 +114,23 @@ def test_very_high_hydrophobic_penalised_more(self): class TestChargeDensityPenalty: def test_charge_density_at_threshold_no_penalty(self): - # |cd| exactly 0.55 → no penalty (condition is > 0.55) + # cd exactly 0.55 → no penalty (condition is > 0.55; abs() removed) assert safety_score(_feat(charge_density=0.55)) == 1.0 - def test_negative_charge_density_uses_abs(self): - # charge_density = -0.60 → |cd| = 0.60 > 0.55 → penalty + def test_negative_charge_density_no_hemolysis_penalty(self): + # Negatively charged peptides are electrostatically repelled from membranes + # and pose no hemolysis risk. Only positive charge_density > 0.55 should trigger + # the penalty. The old abs() behavior was wrong. score_neg = safety_score(_feat(charge_density=-0.60)) score_pos = safety_score(_feat(charge_density=0.60)) - assert score_neg == score_pos + assert score_neg == 1.0, ( + f"Negatively charged peptide (cd=-0.60) should have no hemolysis penalty, " + f"got safety_score={score_neg}" + ) + assert score_pos < 1.0, ( + f"Highly positively charged peptide (cd=0.60) should have hemolysis penalty, " + f"got safety_score={score_pos}" + ) def test_charge_density_0_70_penalty_correct(self): # cd = 0.70 → risk = (0.70 - 0.55) * 1.2 = 0.18 → score = 0.82