tooluniverse-protein-sae-variant-interpretation

What this skill does

# Protein SAE Variant Interpretation

Interpret a single missense variant by comparing reference vs mutant Sparse Autoencoder (SAE) feature activations from the ESMC-6B protein language model. SAE features are interpretable latent dimensions of the model's hidden state — many activate on biologically meaningful patterns (active sites, ligand-binding pockets, PTM sequons, structural motifs).

---

## When to use this skill

Apply when users:
- Ask "why is variant X (a missense) loss-of-function?" and need a mechanistic answer beyond a pathogenicity score
- Have an AlphaMissense / ClinVar "damaging" variant and want to know **which functional feature breaks** (catalytic? binding? PTM site? structural?)
- Want to compare ref vs mutant protein representation at a specific residue
- Are interpreting why a structurally subtle change (single AA) has a big functional impact

**Not for** (use other skills instead):
- ACMG pathogenicity classification → `tooluniverse-variant-interpretation`
- Regulatory / non-coding variants → `tooluniverse-variant-to-mechanism`
- Variant-to-disease association without mechanism → `tooluniverse-gene-disease-association`
- Cancer-specific variant interpretation → `tooluniverse-cancer-variant-interpretation`

---

## Required inputs

| Input | Format | Example |
|---|---|---|
| Protein identifier | UniProt accession or HGNC gene symbol | `P04637` or `TP53` |
| Variant | Single-letter code: `{ref_aa}{position_1idx}{alt_aa}` | `R175H` |

Optional:
- Window radius (default 8): residues around the mutation to analyze
- Reference protein sequence (skip the UniProt lookup if already known)

---

## Prerequisites

- **ESM_API_KEY** env var with a valid EvolutionaryScale Forge token (https://forge.evolutionaryscale.ai)
- **esm package with SAE support**:
  ```
  pip install 'esm @ git+https://github.com/evolutionaryscale/esm@ee891c52'
  ```
  The PyPI release of `esm` does NOT yet include SAEConfig. Install from the upstream feature branch.

**License note**: SAE outputs from Forge are governed by the [Cambrian Inference Clickthrough License](https://www.evolutionaryscale.ai/policies/cambrian-inference-clickthrough-license-agreement) — non-commercial / academic research use only unless a separate commercial agreement applies.

---

## Workflow (5 steps)

### Step 1: Resolve gene → UniProt accession (if needed)

If the user gave a gene symbol, resolve to a reviewed human accession:

```python
UniProt_search(
    query="gene:TP53 AND organism_id:9606 AND reviewed:true",
    fields=["accession", "gene_names", "protein_name"],
)
# → returns accession P04637 as the canonical reviewed human TP53
```

### Step 2: Fetch the canonical reference sequence

```python
UniProt_get_sequence_by_accession(accession="P04637")
# → returns the canonical isoform 1 sequence (393 AA for TP53)
```

### Step 3: Validate the reference residue + build mutant sequence

Parse the variant string (e.g. `R175H`):
- `ref_aa = "R"`, `position = 175` (1-indexed), `alt_aa = "H"`
- Verify `ref_sequence[174] == "R"` (Python is 0-indexed, position is 1-indexed)
- Build mutant: `mutant_sequence = ref_sequence[:174] + "H" + ref_sequence[175:]`

If the reference residue does NOT match, return an explicit error — do not silently mutate the wrong position.

### Quick path (recommended for variant analysis): composite tool

For the standard variant-interpretation use case, use one of two composite tools depending on how much you need:

**Fullest one-call** — disruption + per-feature biological category labels + mechanism summary:
```python
ESM_explain_variant_mechanism(
    sequence=ref_sequence,
    position=175, ref_aa="R", alt_aa="H",
    window=8,
    top_k_features=5,       # describe top 5 lost + top 5 gained
)
# data["mechanism_summary"] e.g. "Disrupted feature categories (lost): catalytic=2, ligand-binding=1"
# data["lost_feature_categories"] / ["gained_feature_categories"] — category counts
# data["top_features_lost"] / ["top_features_gained"] — per-feature delta + category + confidence
```
This is the right default for variant-mechanism reports — saves you from chaining `ESM_score_variant_sae_disruption` + N `ESM_describe_sae_feature` calls. Set `include_descriptions=false` to skip labeling (2 Forge calls only) when you just need the deltas.

**Raw delta only** (no category labels, no describe calls — faster):
```python
ESM_score_variant_sae_disruption(
    sequence=ref_sequence,
    position=175, ref_aa="R", alt_aa="H",
    window=8, top_k_features=10,
)
# → returns top_features_lost + top_features_gained ranked by |delta|
#   plus ref / mut activation sums per feature
```

If ref_aa doesn't match the sequence at the given position, both tools return a clear error (you supplied the wrong isoform / mis-labeled the variant). The longer path below is for inspecting raw per-residue features.

**Multiple variants at once** (e.g. saturation at residue 175 — all 19 alternates):
```python
alts = "ACDEFGHIKLMNPQRSTVWY".replace("R", "")
variants = [{"position": 175, "ref_aa": "R", "alt_aa": a} for a in alts]
ESM_score_variant_sae_batch(sequence=ref_sequence, variants=variants, top_k_features=5)
# 1 + 19 = 20 Forge calls, not 38
```

### Step 4 (long path): Get SAE features for reference and mutant

```python
ref_features = ESM_get_sae_features(
    sequence=ref_sequence,
    position=175,          # 1-indexed mutation position
    window=8,              # +/- 8 residues = 17-residue window
    top_k_per_residue=64,  # full sparsity (k=64 is the SAE's actual k)
)

mut_features = ESM_get_sae_features(
    sequence=mutant_sequence,
    position=175,
    window=8,
    top_k_per_residue=64,
)
```

Each call returns a list of `{residue_idx_1based, active_features: [{feature_id, activation}]}` for residues in the window. Typical latency: ~1-3 seconds per call (so ~2-6s total). Forge cost: 2 credits (1 per call).

### Step 5: Compute per-feature activation deltas

Aggregate gained / lost features over the window:

```python
# Build per-feature activation arrays across the window
def feature_to_window_sum(features_response):
    sums = {}  # feature_id -> sum of activations across all residues in window
    for residue in features_response["data"]["activations"]:
        for f in residue["active_features"]:
            sums[f["feature_id"]] = sums.get(f["feature_id"], 0.0) + f["activation"]
    return sums

ref_sums = feature_to_window_sum(ref_features)
mut_sums = feature_to_window_sum(mut_features)

# Delta = mut - ref. Positive = gained on mutation. Negative = lost.
all_features = set(ref_sums) | set(mut_sums)
deltas = {
    f: mut_sums.get(f, 0.0) - ref_sums.get(f, 0.0)
    for f in all_features
}

top_lost = sorted(deltas.items(), key=lambda x: x[1])[:10]       # most negative
top_gained = sorted(deltas.items(), key=lambda x: -x[1])[:5]     # most positive
```

---

## Interpretation table

The 16,384 SAE features have been categorized (by UniRef90 activation patterns) into ~6 biological types. Interpret features by their dominant category:

| Category | What it means biologically | Variant types it often catches |
|---|---|---|
| **Catalytic function** | Activates on residues at or near enzyme active sites | Variants that disrupt enzymatic activity (kinases, hydrolases, transferases) |
| **Ligand-binding site** | Activates on residues that contact small-molecule / ion / nucleotide ligands | Variants disrupting drug binding, ATP/GTP binding, metal coordination, DNA/RNA binding |
| **Post-translational modification (PTM)** | Activates on phospho-sites, glycosylation sequons, ubiquitin sites, acetylation sites | Variants disrupting phospho-regulation, N-glycosylation, ubiquitin-mediated degradation |
| **Domain / motif** | Activates on classic structural domains (Zn finger, leucine zipper, EF-hand, etc.) | Variants disrupting tertiary fold within a domain |
| **Structural stability** | Activates on residues critical to local fold | Variants destabilizing the protein → folding LoF |
|