diff --git a/Makefile b/Makefile index 09966abd..a44b5eeb 100644 --- a/Makefile +++ b/Makefile @@ -20,7 +20,7 @@ lint: $(RUFF) check src tests scripts ci: lint test - @echo "CI passed: lint + 1088-test suite green" + @echo "CI passed: lint + 1124-test suite green" bench-leakage: PYTHONPATH=src $(PYTHON) -m openamp_foundry.cli bench leakage \ diff --git a/src/openamp_foundry/features/physchem.py b/src/openamp_foundry/features/physchem.py index ef2b9482..88e33302 100644 --- a/src/openamp_foundry/features/physchem.py +++ b/src/openamp_foundry/features/physchem.py @@ -13,6 +13,11 @@ # Trypsin cleaves after K or R; chymotrypsin after F, W, Y (interior sites only) TRYPSIN_SITES = set("KR") CHYMOTRYPSIN_SITES = set("FWY") # identical to AROMATIC — kept separate for semantic clarity +# Human neutrophil elastase (HNE) primary P1 substrates: Ala > Val > Ser +# Reference: Bieth (1986) Bull Eur Physiopathol Respir; Doherty et al. (1991) Biochemistry +ELASTASE_SITES = {"A", "V", "S"} +# 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. @@ -181,6 +186,52 @@ def selectivity_proxy(charge: float, gravy: float) -> float: return round(0.6 * charge_sel + 0.4 * gravy_sel, 4) +def aggregation_propensity(sequence: str) -> float: + """Heuristic aggregation propensity for short synthetic peptides [0, 1]. + + Two-component model: + 1. Hydrophobic run length: consecutive residues from {V,I,L,M,F,W} ≥ 4 drive + self-association during SPPS, solubilisation, and assay buffers. + Risk onset at run ≥ 4 (same threshold as QC HYDROPHOBIC_RUN_RE flag; same + residue set AGG_HYDROPHOBIC). Ramp: 0.0 at run=4 → 1.0 at run ≥ 8. + 2. Beta-branched residue density (Val, Ile, Thr) > 20% promotes β-strand over + α-helix and intermolecular β-sheet aggregation in concentrated solutions. + + Limitation: Ala-rich sequences score 0.0 despite documented on-resin aggregation + issues (Sarin et al. 1984; Tam et al. 1988). Ala aggregation arises from apolar + solvation collapse rather than the hydrophobic-run / β-branched mechanisms here. + + References: + - Quittot et al. (2017) Protein Sci 26:720-735 (hydrophobic run aggregation) + - Wurth et al. (2006) J Mol Biol 355:524-536 (Val/Ile beta-strand aggregation) + + Returns 0.0 (no risk) to 1.0 (severe aggregation risk). + """ + if not sequence: + return 0.0 + n = len(sequence) + + # Component 1: longest hydrophobic run across full sequence (mirrors QC regex) + max_run = 0 + current_run = 0 + for aa in sequence: + if aa in AGG_HYDROPHOBIC: + current_run += 1 + max_run = max(max_run, current_run) + else: + current_run = 0 + # Ramp: 0.0 at run < 4, 0.20 at run=4, …, 1.0 at run ≥ 8 + run_risk = min(1.0, max(0.0, (max_run - 3) / 5.0)) + + # Component 2: beta-branched residue density (Val, Ile, Thr) + beta_branched = {"V", "I", "T"} + beta_density = sum(1 for aa in sequence if aa in beta_branched) / n + # Risk onset at >20%, full risk at >50% + beta_risk = min(1.0, max(0.0, (beta_density - 0.20) / 0.30)) if beta_density > 0.20 else 0.0 + + return round(0.7 * run_risk + 0.3 * beta_risk, 4) + + def compute_features(sequence: str) -> dict[str, float | int | dict[str, int]]: counts = Counter(sequence) length = len(sequence) @@ -194,13 +245,16 @@ def compute_features(sequence: str) -> dict[str, float | int | dict[str, int]]: mu_h = 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) # site density: interior cleavage sites per residue (0 = stable, 1 = all residues cleave) trypsin_density = round(n_trypsin / length if length else 0.0, 4) chymotrypsin_density = round(n_chymotrypsin / length if length else 0.0, 4) + elastase_density = round(n_elastase / length if length else 0.0, 4) charge_ph74 = net_charge_at_ph74(sequence) helix_pa = helix_propensity_score(sequence) gravy = gravy_score(sequence) sel_proxy = selectivity_proxy(charge_ph74, gravy) + agg = aggregation_propensity(sequence) return { "length": length, "net_charge_proxy": charge, @@ -218,9 +272,12 @@ def compute_features(sequence: str) -> dict[str, float | int | dict[str, int]]: "boman_index": boman_index(sequence), "gravy": gravy, "selectivity_proxy": sel_proxy, + "aggregation_propensity": agg, "residue_counts": dict(sorted(counts.items())), "trypsin_site_density": trypsin_density, "chymotrypsin_site_density": chymotrypsin_density, + "elastase_site_density": elastase_density, "interior_trypsin_sites": n_trypsin, "interior_chymotrypsin_sites": n_chymotrypsin, + "interior_elastase_sites": n_elastase, } diff --git a/src/openamp_foundry/scoring/stability.py b/src/openamp_foundry/scoring/stability.py index cbe9423f..34379710 100644 --- a/src/openamp_foundry/scoring/stability.py +++ b/src/openamp_foundry/scoring/stability.py @@ -6,37 +6,45 @@ def serum_stability_score(features: dict[str, Any]) -> float: - """Predict relative serum stability from interior protease-site density. + """Predict relative stability against serum/tissue proteases from site density. - Short cationic AMPs are degraded primarily by serum trypsin and chymotrypsin. - Each interior K/R site adds trypsin susceptibility; F/W/Y sites add chymotrypsin - susceptibility. Higher site density → lower predicted stability. + Short cationic AMPs face three major proteolytic threats in biological fluids: + 1. Serum trypsin (K/R sites): dominant degradation route in blood/plasma. + 2. Serum chymotrypsin (F/W/Y sites): secondary serum route. + 3. Neutrophil elastase (A/V/S sites): abundant at infection sites; degrades + helix-forming AMPs with high Ala content significantly faster than serum + trypsin alone predicts. + + Weights reflect relative proteolytic efficiency for short cationic peptides: + trypsin 2.0 > chymotrypsin 1.0 > elastase 0.5 Literature basis: - Hilpert et al. (2006) J Antimicrob Chemother: short cationic peptides with ≥3 interior K/R sites typically have serum t½ < 30 min. - Wade et al. (1990) PNAS: D-amino acid substitution extends t½ by >10×. + - Bieth (1986) Bull Eur Physiopathol Respir; Doherty et al. (1991) Biochemistry: + HNE cleaves Ala > Val > Ser; abundant at infection sites (>1 μM). Returns a value in [0, 1]: - 1.0 = no interior cleavage sites (e.g., pure alanine or hydrophobic peptide) - 0.5 = moderate site density (~0.25 per residue) - 0.0 = very high site density (≥0.5 per residue) + 1.0 = no interior cleavage sites + 0.5 = moderate combined site density (~0.25 per residue) + 0.0 = very high combined site density (≥0.5 per residue) - This score is informational — it does not currently gate the ensemble. - Use it to stratify candidates by predicted serum longevity and to flag - candidates that require stability engineering before clinical translation. + This score is informational — it does not gate the ensemble. + Use it to stratify candidates by predicted proteolytic longevity. """ length = features.get("length", 0) if not length: return 1.0 trypsin_density = features.get("trypsin_site_density", 0.0) - # Both densities are pre-normalised by compute_features(); no division needed here. chymo_density = features.get("chymotrypsin_site_density", 0.0) + # Elastase density: pre-normalised by compute_features(); defaults 0 for backward compat. + elastase_density = features.get("elastase_site_density", 0.0) - # Trypsin is the primary serum protease for cationic peptides; weight 2:1 vs chymotrypsin - combined_density = (2.0 * trypsin_density + 1.0 * chymo_density) / 3.0 + # Weighted sum: trypsin 2 : chymotrypsin 1 : elastase 0.5 (denominator = sum of weights) + combined_density = (2.0 * trypsin_density + 1.0 * chymo_density + 0.5 * elastase_density) / 3.5 - # At combined_density = 0.5 (every other residue is a cleavage site), score → 0 + # At combined_density ≈ 0.5 (every other residue is a cleavage site), score → 0 score = 1.0 - clamp01(combined_density / 0.5) return round(score, 4) diff --git a/src/openamp_foundry/scoring/synthesis.py b/src/openamp_foundry/scoring/synthesis.py index ee8e6bab..26c269c6 100644 --- a/src/openamp_foundry/scoring/synthesis.py +++ b/src/openamp_foundry/scoring/synthesis.py @@ -4,11 +4,23 @@ def synthesis_feasibility_score(features: dict, valid_sequence: bool = True) -> float: + """Estimate solid-phase synthesis difficulty [0, 1]; higher = easier to synthesise. + + Penalty sources: + - Length > 30: longer chains accumulate deletion errors and solubility issues. + - Length < 8: too short to be reliably purified / characterised. + - Longest single-AA repeat run ≥ 5: coupling efficiency drops on homo-repeat stretches. + - Cysteine fraction > 20%: disulphide scrambling and side-chain protection cost. + - Aggregation propensity > 0: interior hydrophobic runs (VILMFW ≥ 4) and high + beta-branched density (Val/Ile/Thr) cause on-resin aggregation and poor solubility. + References: Quittot et al. (2017) Protein Sci; Wurth et al. (2006) J Mol Biol. + """ if not valid_sequence: return 0.0 length = features["length"] repeat_run = features["longest_repeat_run"] cys = features["cysteine_fraction"] + agg = features.get("aggregation_propensity", 0.0) score = 1.0 if length > 30: @@ -19,4 +31,6 @@ def synthesis_feasibility_score(features: dict, valid_sequence: bool = True) -> score -= 0.10 if cys > 0.20: score -= 0.15 + if agg > 0.0: + score -= min(agg * 0.25, 0.20) return round(clamp01(score), 4) diff --git a/tests/test_aggregation_propensity.py b/tests/test_aggregation_propensity.py new file mode 100644 index 00000000..3624bcc7 --- /dev/null +++ b/tests/test_aggregation_propensity.py @@ -0,0 +1,133 @@ +"""Tests for aggregation propensity scoring and its synthesis score integration.""" +from __future__ import annotations + +import pytest + +from openamp_foundry.features.physchem import aggregation_propensity, compute_features +from openamp_foundry.scoring.synthesis import synthesis_feasibility_score + + +class TestAggregationPropensityUnit: + def test_empty_sequence_returns_zero(self): + assert aggregation_propensity("") == 0.0 + + def test_polar_sequence_no_aggregation(self): + # All charged/polar residues — no hydrophobic run + assert aggregation_propensity("KKKKKK") == 0.0 + assert aggregation_propensity("EEEEEE") == 0.0 + assert aggregation_propensity("RRRRRRR") == 0.0 + + def test_single_hydrophobic_no_run(self): + # Only one hydrophobic in a polar context — no run ≥ 4 + assert aggregation_propensity("KKVKK") == 0.0 + + def test_short_run_below_threshold(self): + # KKVLLK: run of 3 interior hydrophobics (V,L,L) — below ≥4 threshold → run_risk=0 + # beta_density: only V = 1/6 ≈ 0.167 < 0.20 → beta_risk=0 + assert aggregation_propensity("KKVLLK") == 0.0 + + def test_run_of_4_triggers_run_risk(self): + # "KKVILL": full-sequence run of 4 (VILL at pos 2-5) → run_risk = (4-3)/5 = 0.20 + # Combined with beta component → total > 0.14 (run_risk alone = 0.7 * 0.20 = 0.14) + score = aggregation_propensity("KKVILL") + assert score > 0.14 + + def test_run_of_4_boundary_exact_run_component(self): + # "KVLLLK": full-sequence run exactly 4 (VLLL) → run_risk = (4-3)/5 = 0.20 + # beta_density: V(1)/6 = 0.167 < 0.20 → beta_risk = 0 + # Total = 0.7 * 0.20 = 0.14 + score = aggregation_propensity("KVLLLK") + assert score == pytest.approx(0.14, abs=0.01) + + def test_run_of_8_saturates_at_max_run_risk(self): + # "KVLLLLLLLK": run of 8 (VLLLLLLL) → run_risk = (8-3)/5 = 1.0 (saturated) + # beta_density: V(1)/10 = 0.10 < 0.20 → beta_risk = 0 + # Total = 0.7 * 1.0 = 0.70 + score = aggregation_propensity("KVLLLLLLLK") + assert score == pytest.approx(0.70, abs=0.01) + + def test_long_hydrophobic_run_high_risk(self): + # Run of 8 interior hydrophobics → max run_risk (capped at 1.0) + score = aggregation_propensity("KVIIIIIIIK") + assert score >= 0.6 + + def test_beta_branched_density_alone(self): + # High Val/Ile/Thr density without a long run + # VIVTVITV: V,I,T at 8/8 positions but no consecutive run of 4 agg residues (T not in AGG_HYDROPHOBIC) + # Actually VIT: V and I are in AGG_HYDROPHOBIC, T is not → runs are broken by T + score = aggregation_propensity("VIVTVITV") + # beta_density = 8/8 = 1.0 → beta_risk = max of (1.0-0.20)/0.30 capped at 1 = 1.0 + # run_risk = 0 (runs of V,I broken by T: max run = 2) + # aggregation = 0.7*0 + 0.3*1.0 = 0.3 + assert score == pytest.approx(0.3, abs=0.05) + + def test_pure_ile_long_run_max_risk(self): + # IIIIIIIIIII: all interior = run of 9 + # run_risk = min(1.0, (9-3)/5) = 1.0; beta_density = 1.0 → beta_risk = 1.0 + score = aggregation_propensity("IIIIIIIIIII") + assert score == pytest.approx(1.0, abs=0.01) + + def test_returns_unit_interval(self): + for seq in ["AAAAAAA", "KWKLFKKIGAVLKVL", "GIGKFLHSAKKFGKAFVGEIMNS", + "FLPLIGRVLSGIL", "IIIIIIII", "KKKKKKK"]: + score = aggregation_propensity(seq) + assert 0.0 <= score <= 1.0, f"score={score} out of [0,1] for {seq}" + + def test_compute_features_includes_key(self): + feats = compute_features("KWKLFKKIGAVLKVL") + assert "aggregation_propensity" in feats + + def test_temporin_like_seed004_no_run_only_beta(self): + # FLPLIGRVLSGIL: longest interior hydrophobic run = 2 (FL or VL) — below ≥4 threshold + # beta_density: I(4),V(7),I(11) = 3/13 = 0.231 → small beta_risk component + score = aggregation_propensity("FLPLIGRVLSGIL") + # run_risk = 0; small beta_risk → total < 0.1 + assert score < 0.1 + # And it's genuinely > 0 because of the beta-branched component + assert score > 0.0 + + +class TestSynthesisFeasibilityWithAggregation: + def test_high_aggregation_penalises_synthesis(self): + feats_high = compute_features("IIIIIIIIIIII") # long Ile run → agg=1.0 + feats_low = compute_features("KWKLFKKIGAVLKVL") # typical AMP, low agg + assert synthesis_feasibility_score(feats_high) < synthesis_feasibility_score(feats_low) + + def test_zero_aggregation_no_penalty(self): + # Synthesis score should only be affected by length/cys/repeat, not aggregation + feats_no_agg = { + "length": 8, "longest_repeat_run": 8, "cysteine_fraction": 0.0, + "aggregation_propensity": 0.0, + } + score = synthesis_feasibility_score(feats_no_agg) + assert score == pytest.approx(0.90, abs=0.01) # only repeat_run≥5 penalty applies + + def test_aggregation_penalty_capped_at_0_20(self): + # Max agg (1.0) → penalty = min(1.0 * 0.25, 0.20) = 0.20 + feats_max_agg = { + "length": 10, "longest_repeat_run": 1, "cysteine_fraction": 0.0, + "aggregation_propensity": 1.0, + } + score = synthesis_feasibility_score(feats_max_agg) + assert score == pytest.approx(0.80, abs=0.01) + + def test_mild_aggregation_modest_penalty(self): + feats_mild = { + "length": 10, "longest_repeat_run": 1, "cysteine_fraction": 0.0, + "aggregation_propensity": 0.4, + } + score = synthesis_feasibility_score(feats_mild) + # penalty = 0.4 * 0.25 = 0.10 + assert score == pytest.approx(0.90, abs=0.01) + + def test_synthesis_score_unit_interval(self): + for seq in ["IIIIIIIIIII", "KWKLFKKIGAVLKVL", "GIGKFLHSAKKFGKAFVGEIMNS", "KKKK"]: + feats = compute_features(seq) + s = synthesis_feasibility_score(feats) + assert 0.0 <= s <= 1.0, f"score={s} out of [0,1] for {seq}" + + def test_backward_compat_without_aggregation_key(self): + # Caller passes features dict without aggregation_propensity key (pre-this-PR dict) + feats_old = {"length": 15, "longest_repeat_run": 1, "cysteine_fraction": 0.0} + score = synthesis_feasibility_score(feats_old) + assert score == 1.0 diff --git a/tests/test_elastase_stability.py b/tests/test_elastase_stability.py new file mode 100644 index 00000000..8e13300f --- /dev/null +++ b/tests/test_elastase_stability.py @@ -0,0 +1,103 @@ +"""Tests for elastase-extended serum stability scoring.""" +from __future__ import annotations + +from openamp_foundry.features.physchem import compute_features, interior_protease_sites, ELASTASE_SITES +from openamp_foundry.scoring.stability import serum_stability_score + + +class TestElastaseSites: + def test_ala_is_elastase_site(self): + assert "A" in ELASTASE_SITES + + def test_val_is_elastase_site(self): + assert "V" in ELASTASE_SITES + + def test_ser_is_elastase_site(self): + assert "S" in ELASTASE_SITES + + def test_lys_not_elastase_site(self): + assert "K" not in ELASTASE_SITES + + def test_phe_not_elastase_site(self): + assert "F" not in ELASTASE_SITES + + def test_pure_ala_interior_sites(self): + # AAAAAAAAA (9 aa): 8 interior sites (excluding C-terminal A) + n = interior_protease_sites("AAAAAAAAA", ELASTASE_SITES) + assert n == 8 + + def test_no_elastase_site_sequence(self): + # KKKFWYKKKR — no A/V/S sites + n = interior_protease_sites("KKKFWYKKKR", ELASTASE_SITES) + assert n == 0 + + def test_cterminal_excluded(self): + # KA: only 1 residue before C-terminal; A is at position 0 (interior), excluded by C-term rule + # sequence "KA": interior_protease_sites excludes C-terminal (index 1=A) + # so only K (index 0) is checked: K not in ELASTASE_SITES → 0 + n = interior_protease_sites("KA", ELASTASE_SITES) + assert n == 0 + + def test_compute_features_has_elastase_keys(self): + feats = compute_features("KWKLFKKIGAVLKVL") + assert "elastase_site_density" in feats + assert "interior_elastase_sites" in feats + + def test_elastase_density_is_unit_interval(self): + for seq in ["AAAAAA", "KWKLFK", "GIGKFLHSAKKFGKAFVGEIMNS", "FLPLIGRVLSGIL"]: + feats = compute_features(seq) + d = feats["elastase_site_density"] + assert 0.0 <= d <= 1.0, f"elastase_density={d} out of [0,1] for {seq}" + + +class TestSerumStabilityWithElastase: + def test_high_ala_sequence_lower_stability_than_lys_rich(self): + # Poly-Ala has many elastase sites; high-Lys has many trypsin sites + # Both should get somewhat reduced stability, elastase modestly penalises Ala + ala_seq = "AAAAAAAAAAAAA" + ala_feats = compute_features(ala_seq) + ala_score = serum_stability_score(ala_feats) + # Pure Ala has NO trypsin or chymotrypsin sites → old score was 1.0 + # Now elastase (A sites) gives it a modest penalty + assert ala_score < 1.0, "Poly-Ala should have reduced stability due to elastase sites" + + def test_no_cleavage_site_sequence_perfect_stability(self): + # Poly-Ile has no trypsin, chymotrypsin, OR elastase sites + feats = compute_features("IIIIIIIIIII") + score = serum_stability_score(feats) + assert score == 1.0 + + def test_old_trypsin_dominant_still_holds(self): + # High Lys (trypsin sites) should still give the lowest stability + high_lys = compute_features("KKKKKKKKKKKK") + high_ala = compute_features("AAAAAAAAAAAA") + assert serum_stability_score(high_lys) < serum_stability_score(high_ala) + + def test_elastase_feature_absent_backward_compat(self): + # If a caller passes a features dict without elastase_site_density (pre-PR features), + # serum_stability_score should still work (default=0.0 for elastase term) + feats = {"length": 15, "trypsin_site_density": 0.1, "chymotrypsin_site_density": 0.0} + score = serum_stability_score(feats) + assert 0.0 <= score <= 1.0 + + def test_elastase_adds_modest_penalty_to_typical_amp(self): + # Magainin-2 analog: GIGKFLHSAKKFGKAFVGEIMNS + # Should have lower stability than a version with all Ala→Ile + magainin = compute_features("GIGKFLHSAKKFGKAFVGEIMNS") + # Replace all A/V/S with I (no elastase sites) + no_elastase = compute_features("GIGKFLHIKKFGKIFVGEIMNS".replace("S", "I").replace("A", "I").replace("V", "I")) + assert serum_stability_score(magainin) <= serum_stability_score(no_elastase) + + def test_score_is_unit_interval(self): + for seq in ["AAAAAA", "KWKLFK", "GIGKFLHSAKKFGKAFVGEIMNS", "FLPLIGRVLSGIL", "RRRRRR"]: + feats = compute_features(seq) + s = serum_stability_score(feats) + assert 0.0 <= s <= 1.0, f"score={s} out of [0,1] for {seq}" + + def test_denominator_normalisation_prevents_unbounded_penalty(self): + # Worst-case: sequence with all K/R (trypsin), all F/Y/W (chymotrypsin), all A/V/S (elastase) + # is impossible (same residues can't be all three), but combined density ≤ 1.0 + feats = {"length": 10, "trypsin_site_density": 1.0, + "chymotrypsin_site_density": 1.0, "elastase_site_density": 1.0} + score = serum_stability_score(feats) + assert score == 0.0