tooluniverse-residue-functional-mechanism-interpretation

What this skill does

# Residue functional mechanism interpretation

The core user question: **"Why are these residues functionally critical?"**

Answering needs more than one source: a residue in the ligand pocket means
something different from one in a catalytic triad, an interface, the
hydrophobic core, or a PTM sequon. This skill synthesizes evidence from
multiple TU tools to call a mechanism for each residue (or cluster of
adjacent residues).

The residues can come from **any source** — this skill is agnostic to where
they came from:

| Residue source | Typical call pattern |
|---|---|
| **DMS map hotspots** | Use Step 1 (optional) to detect top-K by max effect, then continue |
| **ClinVar recurrent variants** | Pull recurrent positions from ClinVar, pass directly as `user_provided_positions` |
| **Literature hot regions** | Paste positions from a paper's Fig 1, pass directly |
| **Evolutionarily conserved residues** | Filter by conservation score, pass top-N positions |
| **Druggable site residues** | From a binding-site predictor, pass directly |
| **Clinician's question** | "Why does mutation at R175 keep showing up in tumors?" — pass `[175]` |

---

## When to use this skill

- You have a list of residues (from anywhere) and need to explain *why* they
  matter biologically
- You're writing the methods/discussion section of a paper and need
  mechanistic claims per residue or cluster
- You're comparing residue sets across orthologs and need a per-residue
  category to align by

**Not for**:
- Validating a *predictor* against DMS — use
  `tooluniverse-variant-predictor-dms-validation`
- Single-variant **SAE feature decomposition** when you want to see which
  ESMC features changed — use `tooluniverse-protein-sae-variant-interpretation`
- Single-variant **LoF mechanism synthesis** for one variant in isolation
  — use `tooluniverse-protein-lof-mechanism`

The single-variant skills focus on one mutation's signature; this skill
focuses on *which residues matter and why*.

---

## Required inputs (choose one entry path)

**Path A — Residues supplied directly** (covers ClinVar / literature /
custom positions):

| Input | Notes |
|---|---|
| `user_provided_positions: List[int]` | 1-based canonical residue positions |
| Protein metadata | UniProt accession + PDB ID + chain |
| (optional) DMS matrix | If supplied, used to enrich each cluster with effect-size context |
| (optional) SAE tensor | If supplied, gives SAE feature evidence as a 4th layer |

**Path B — Detect hotspots from a DMS map** (original use case):

| Input | Source | Notes |
|---|---|---|
| DMS effect matrix `(20, n_positions)` | `MaveDB_get_effect_matrix` | NaN for unmeasured |
| `disruptive_tail` | DMS retrieval metadata | `"top"` or `"bottom"` |
| Protein metadata | UniProt accession + PDB ID + chain | for the multi-evidence lookups |
| (optional) SAE evidence | `ESM_get_region_sae_features` for one contiguous cluster (1 Forge call), OR a precomputed full DMS SAE tensor from `ESM_get_sae_features` per mutant | the SAE evidence layer; see Step 4 for which path |

In Path B the skill runs Step 1 to detect hotspots; in Path A it skips Step
1 entirely and goes straight to Step 2 (gather evidence).

---

## Workflow

### Step 0 (MANDATORY if user names specific positions): Premise check

If the user says "explain why residue/cluster X is a hotspot", do NOT take
that as a given. Verify it's actually a hotspot in THIS DMS first — users
import biological knowledge from other contexts that may not match what the
specific assay measured.

```python
# Per-position disruption magnitude (same formula as Step 1)
if disruptive_tail == "top":
    dms_per_pos = np.nanmax(dms_matrix, axis=0)
elif disruptive_tail == "bottom":
    dms_per_pos = -np.nanmin(dms_matrix, axis=0)

# Where do the user-named positions actually rank?
ranks = (-dms_per_pos).argsort().argsort()  # 0 = highest
for user_pos in user_named_positions:
    rank = int(ranks[pos_index[user_pos]])
    pct = 100 * (1 - rank / len(dms_per_pos))
    print(f"  pos {user_pos}: rank {rank+1}/{len(dms_per_pos)}, "
          f"top {pct:.0f}% by max effect")
```

**Decision rule:**
- If named position is in top 25% by max effect → premise confirmed, proceed.
- If named position is in top 50% but not 25% → premise weakly supported;
  proceed but note the rank in your report.
- **If named position is below top 50% → REPORT THIS MISMATCH TRANSPARENTLY**
  at the top of your answer before continuing with the mechanism analysis.

Concrete example: the user asks "why is KRAS G12/G13 a folding hotspot in
this DMS?" The data say G12 ranks 105/187 (top 56%) and G13 ranks 124/187
(top 66%) by max ΔΔG — they are NOT folding hotspots in this AbundancePCA
assay (even though they ARE famous oncogenic positions). The skill's job
is to surface that contradiction up front, then proceed with a mechanism
analysis for those residues (their oncogenic effect is via GTPase
abolishment, not fold disruption — and that's a genuinely useful answer to
the user's actual scientific question, just not the one they literally asked).

### Step 1 (Path B only): Detect hotspots from the DMS matrix

**Skip this step entirely if the user already provided positions (Path A).**
In Path A the residue list IS the input; jump straight to clustering at the
bottom of this section.

```python
import numpy as np

if user_provided_positions:
    # Path A — residues from any source (ClinVar / literature / custom)
    positions_to_analyze = sorted(set(user_provided_positions))
else:
    # Path B — detect from DMS matrix
    if disruptive_tail == "top":
        dms_per_pos = np.nanmax(dms_matrix, axis=0)   # max destabilization at any allele
    elif disruptive_tail == "bottom":
        dms_per_pos = -np.nanmin(dms_matrix, axis=0)  # flip for low-is-bad assays

    K = 20
    positions_to_analyze = sorted(np.argsort(-dms_per_pos)[:K].tolist())

# Chain adjacent positions (gap ≤ 2) into clusters — shared by both paths
clusters = []
current = [positions_to_analyze[0]]
for p in positions_to_analyze[1:]:
    if p - current[-1] <= 2:
        current.append(p)
    else:
        clusters.append(current)
        current = [p]
clusters.append(current)
print(f"{len(clusters)} cluster(s): {clusters}")
```

A "cluster" of one is allowed — it just gets less statistical power in the
permutation test, but multi-evidence interpretation still works.

**Path A workflow contracts** (what's available vs not):
| Evidence layer | Path A (user residues) | Path B (DMS hotspots) |
|---|---|---|
| Structural (Step 2) | ✓ always | ✓ always |
| UniProt features (Step 3) | ✓ always | ✓ always |
| SAE per-feature labels (Step 4, descriptive) | ✓ if SAE tensor supplied | ✓ if SAE tensor supplied |
| SAE permutation test (Step 4-alt) | ✗ — needs DMS-derived `max_drop` baseline; not meaningful for residues with no DMS context | ✓ if SAE tensor supplied |
| DMS effect-size context | ✓ if DMS matrix supplied (enrichment only) | ✓ always |
| Mechanism synthesis (Step 5) | ✓ always | ✓ always |

### Step 2: Gather structural evidence per cluster

Annotate the protein structure once, then read fields for each cluster's
positions:

```python
# One-shot structural annotation (cached for the rest of the skill)
struct = Structure_annotate_per_residue(
    pdb_id="6VJJ",                 # pick a structure with the relevant complex
    target_chain="A",
    partner_chains=["B"],          # if there's a binding partner
    ligand_resnames=["GNP", "MG"], # if there's a relevant ligand
    distance_cutoff=5.0,
    include_secondary_structure=True,
)
by_pos = {r["position"]: r for r in struct["data"]["annotations"]}

for cluster in clusters:
    structural_summary = {
        "interface_count": sum(1 for p in cluster if by_pos.get(p, {}).get("region") in ("interface", "both")),
        "ligand_pocket_count": sum(1 for p in cluster if by_pos.get(p, {}).get("region") in ("ligand", "both")),
        "core_count": sum(1 for p in cluster if by_pos.get(p, {}).get("is