tooluniverse-protein-lof-mechanism

What this skill does

# Protein LoF Mechanism Synthesis

For a single missense variant, integrate 5 independent computational signals to propose a specific loss-of-function mechanism. Each signal answers a different question:

| Signal | Tool | Answers |
|---|---|---|
| Pathogenicity | `AlphaMissense_get_variant_score` | "Is this variant damaging?" |
| Structural context | `alphafold_get_prediction` | "Is the mutation in a folded vs disordered region?" |
| Sequence likelihood | `ESM_score_sequence` | "Is the substitution evolutionarily plausible?" |
| Feature disruption | `ESM_explain_variant_mechanism` (or `ESM_score_variant_sae_disruption` + `ESM_describe_sae_feature` for raw control) | "Which biological feature breaks?" |
| Stability | `DynaMut2_predict_stability` | "Does the protein still fold correctly?" |

---

## When to use this skill

Apply for **missense (coding) variants** where:
- You already have evidence the variant is damaging (or want to find out)
- You need to know **WHY** it's damaging in molecular terms
- Downstream work depends on the mechanism: drug rescue strategies need to know what's broken, gene-therapy decisions need to distinguish "fix the protein" vs "replace the protein", clinical reporting wants a mechanism narrative

**Not for** (use other skills instead):
- Non-coding / regulatory variants → `tooluniverse-variant-to-mechanism`
- ACMG pathogenicity classification → `tooluniverse-variant-interpretation`
- Just the SAE feature disruption (without full synthesis) → `tooluniverse-protein-sae-variant-interpretation`
- Cancer-specific drivers → `tooluniverse-cancer-variant-interpretation`

---

## Required inputs

| Input | Format | Example |
|---|---|---|
| Variant ID | `{accession}_{ref_aa}{position}{alt_aa}` | `P04637_R175H` |
| Or: UniProt accession + variant string | accession + ref/pos/alt | `P04637`, `R175H` |
| Or: gene symbol + variant | HGNC symbol + ref/pos/alt | `TP53`, `R175H` |

Parse the variant string:
```python
import re
m = re.match(r"([A-Z])(\d+)([A-Z])", variant_str)
ref_aa, position, alt_aa = m.group(1), int(m.group(2)), m.group(3)
```

---

## Prerequisites

- **ESM_API_KEY** env var for SAE signals (https://forge.evolutionaryscale.ai)
- **esm package with SAE support**: `pip install 'esm @ git+https://github.com/evolutionaryscale/esm@ee891c52'`
- Internet access for AlphaMissense, AlphaFold, UniProt, DynaMut2 (all hosted endpoints, no API key required for these)

---

## Workflow

### Step 0: Resolve accession + fetch canonical sequence

If user gave a gene symbol:
```python
UniProt_search(
    query="gene:TP53 AND organism_id:9606 AND reviewed:true",
    fields=["accession"],
)
# → P04637
```

Get the canonical sequence:
```python
UniProt_get_sequence_by_accession(accession="P04637")
```

Validate the reference residue:
```python
assert sequence[position - 1] == ref_aa, "Wrong sequence or wrong isoform"
```

### Step 1: AlphaMissense pathogenicity

```python
AlphaMissense_get_variant_score(
    uniprot_id="P04637",
    position=175,
    ref_aa="R",
    alt_aa="H",
)
# → returns score 0..1; ≥0.564 is the "likely pathogenic" threshold per the
#   AlphaMissense paper. Above 0.9 is very confident damaging.
```

If AlphaMissense says benign (≤0.34), the rest of the analysis is exploratory — most benign variants don't have a clear LoF mechanism.

### Step 2: AlphaFold structural context

```python
alphafold_get_prediction(uniprot_id="P04637")
# → returns structure data including per-residue pLDDT
```

Check pLDDT at the mutation position:
- **pLDDT > 70** → well-folded region; the variant is in a structured area; structural / functional disruption is meaningful
- **pLDDT 50-70** → flexible / partially folded; interpretation is ambiguous
- **pLDDT < 50** → disordered; SAE / stability signals may not be reliable

### Step 3: ESMC sequence likelihood

```python
ESM_score_sequence(
    sequence=ref_sequence,
    model="esmc-600m-2024-12",  # or 300m for cheaper
)
# Then for the mutant sequence:
ESM_score_sequence(
    sequence=mutant_sequence,
    model="esmc-600m-2024-12",
)
```

Compute **ΔlogP = mean_logP(mutant) − mean_logP(reference)** at the position window:
- **ΔlogP < −1** → mutation is evolutionarily implausible (strong signal of functional cost)
- **ΔlogP ≈ 0** → mutation is conservative or in a tolerant position
- ΔlogP > 0 is rare and usually noise

### Step 4: SAE feature disruption (the unique signal)

**Recommended — one call** (composite tool, returns disruption + category labels + summary):
```python
ESM_explain_variant_mechanism(
    sequence=ref_sequence,
    position=175,
    ref_aa="R",
    alt_aa="H",
    window=8,
    top_k_features=5,
)
# data["mechanism_summary"] e.g.:
#   "Disrupted feature categories (lost): catalytic=2, ligand-binding=1"
# data["lost_feature_categories"] / data["gained_feature_categories"] give the raw counts
# data["top_features_lost"] / data["top_features_gained"] include per-feature deltas + categories
```

**Lower-level alternative** (use only if you need the raw feature_ids before labeling, e.g. to filter to one category before describe calls):
```python
ESM_score_variant_sae_disruption(
    sequence=ref_sequence, position=175, ref_aa="R", alt_aa="H",
    window=8, top_k_features=10,
)
# Then for each kept feature:
ESM_describe_sae_feature(feature_id=feat["feature_id"])
# → returns category: catalytic | ligand-binding | ptm | domain |
#   motif | structural-stability | secondary-structure |
#   transmembrane | signal-peptide | propeptide | uncategorized
```

**Dominant category among top lost features = the function most likely disrupted.**

### Step 5: DynaMut2 stability (ΔΔG)

DynaMut2 needs a PDB structure. Two options:

**Option A** — use a PDB ID if available for this protein:
```python
# Look up PDB cross-references from UniProt entry
UniProt_get_entry_by_accession(accession="P04637")
# → check `uniProtKBCrossReferences` for entries with database == "PDB"
# → pick a structure that covers the mutation position
```

```python
DynaMut2_predict_stability(
    pdb_id="2FEJ",   # example TP53 DNA-binding domain crystal structure
    chain="A",
    mutation="R175H",
)
# → returns ddG in kcal/mol
```

**Option B** — if no experimental PDB covers the position, use the AlphaFold model (output of Step 2). DynaMut2 accepts AlphaFold PDBs the same way.

### Step 6: Synthesis — decide the LoF mechanism category

Apply the upstream variant_lof_mechanism decision rule:

| Signal pattern | Inferred mechanism |
|---|---|
| ddG > +1 kcal/mol **AND** ΔlogP < 0 | **Structural stability LoF** — mutation destabilizes the fold; protein may misfold / be degraded. Drug rescue strategy: pharmacological chaperones, refolding agents. |
| ddG ≈ 0 (in [-0.5, +1]) **AND** SAE features lost are catalytic | **Direct catalytic LoF** — protein folds normally but the active site is broken. Strategy: substrate analog / cofactor supplementation. |
| ddG ≈ 0 **AND** SAE features lost are ligand-binding | **Binding LoF** — fold preserved, binding pocket disrupted. Strategy: small-molecule restoration. |
| ddG ≈ 0 **AND** SAE features lost are PTM | **PTM LoF** — regulatory site (phospho / glyco / ubiquitin) broken. Mechanism: dysregulation, not direct activity loss. |
| ddG ≈ 0 **AND** SAE features lost are domain / motif | **Interface LoF** — protein-protein interaction surface affected. Strategy: PPI restoration. |
| ddG > 0 + AlphaMissense pathogenic + ΔlogP < 0 but no clear SAE signal | **Generic damaging mutation** — clearly bad but mechanism unclear. Investigate via experimental assay. |

### Step 7: Honest evidence grading

Before reporting, score the synthesis:

| Confidence | Signal requirement |
|---|---|
| **High** | ≥4 signals point the same direction (e.g. AlphaMissense pathogenic + low ΔlogP + ddG > +1 + SAE feature loss) |
| **Medium** | 2-3 signals agree but 1+ are inconclusive |
| **Low** | Signals conflict (e.g. AlphaMissense pathogenic but SAE shows no specific category) — flag for experimental fo