Given a set of residues in a protein, explain WHY they are functionally critical by combining structural context (binding interface, ligand pocket, core, secondary structure), UniProt features (active sites, binding sites, PTM sites, disulfides), optional SAE feature evidence, and optional DMS data. Accepts residues from any source: DMS hotspots (top-K by max effect), ClinVar recurrent variants, literature-reported hot regions, evolutionarily conserved positions, or user-curated lists. Returns a per-cluster mechanism call: catalytic / ligand-binding / interface / structural-core / PTM / regulatory / unknown.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
Given a set of residues in a protein, explain WHY they are functionally critical by combining structural context (binding interface, ligand pocket, core, secondary structure), UniProt features (active sites, binding sites, PTM sites, disulfides), optional SAE feature evidence, and optional DMS data. Accepts residues from any source: DMS hotspots (top-K by max effect), ClinVar recurrent variants, literature-reported hot regions, evolutionarily conserved positions, or user-curated lists. Returns a per-cluster mechanism call: catalytic / ligand-binding / interface / structural-core / PTM / regulatory / unknown.
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.
# 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 = highestfor 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.
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 matrixif disruptive_tail == "top":
dms_per_pos = np.nanmax(dms_matrix, axis=0) # max destabilization at any alleleelif 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:
# 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(1for p in cluster if by_pos.get(p, {}).get("region") in ("interface", "both")),
"ligand_pocket_count": sum(1for p in cluster if by_pos.get(p, {}).get("region") in ("ligand", "both")),
"core_count": sum(1for p in cluster if by_pos.get(p, {}).get("is_core")),
"ss_elements": [by_pos.get(p, {}).get("ss_element") for p in cluster],
}
Step 3: Gather UniProt feature evidence per cluster
# UniProt features per position (active sites, binding sites, PTM sites, etc.)
up = UniProt_get_function_by_accession(accession="P01116")
# Parse up["features"] into a per-position lookup
features_by_pos = {}
for feat in up.get("features", []):
start = int(feat.get("begin", feat.get("position", -1)))
end = int(feat.get("end", start))
for p inrange(start, end + 1):
features_by_pos.setdefault(p, []).append({
"type": feat.get("type"),
"description": feat.get("description", ""),
})
for cluster in clusters:
uniprot_summary = {
"active_site": [p for p in cluster ifany(f["type"] == "Active site"for f in features_by_pos.get(p, []))],
"binding_site": [p for p in cluster ifany(f["type"] == "Binding site"for f in features_by_pos.get(p, []))],
"modified_residue": [p for p in cluster ifany(f["type"] == "Modified residue"for f in features_by_pos.get(p, []))],
"disulfide_bond": [p for p in cluster ifany(f["type"] == "Disulfide bond"for f in features_by_pos.get(p, []))],
"domain": [f.get("description") for p in cluster for f in features_by_pos.get(p, []) if f.get("type") == "Domain"],
}
Two paths depending on whether you already have a full DMS SAE tensor:
Path A โ no precomputed tensor (most cases): use the region tool directly.
If the cluster is a contiguous range (or you can pad to one), ESM_get_region_sae_features aggregates SAE features over the range in a single Forge call:
region = ESM_get_region_sae_features(
sequence=ref_sequence,
start_position=min(cluster_positions),
end_position=max(cluster_positions),
top_k_features=5,
)
top_features = [f["feature_id"] for f in region["data"]["top_features"]]
This is the right default โ 1 Forge call vs 20 ร cluster-size for the DMS-tensor path. For non-contiguous clusters, run once per contiguous sub-range and union the top-K.
Path B โ you already have a precomputed SAE tensor from a DMS sweep (e.g. the variant-predictor-dms-validation pipeline left one on disk): compute drops directly without re-calling Forge.
defcluster_sae_features(sae_tensor, wt_vec, cluster_positions, top_n=5):
"""Returns top SAE features by mean drop at cluster, ready for labeling."""
drops = np.maximum(0.0, wt_vec[None, :, :] - sae_tensor)
max_drop_per_pos = np.nanmax(drops, axis=0) # (n_pos, 16384)
cluster_mean = max_drop_per_pos[cluster_positions].mean(axis=0)
return np.argsort(-cluster_mean)[:top_n].tolist()
top_features = cluster_sae_features(sae_tensor, wt_vec, cluster_cols, top_n=5)
Label each top feature via the SAE feature labeler:
for f in top_features:
label = ESM_describe_sae_feature(feature_id=int(f), n_proteins=5)
print(f" feature {f}: {label['data'].get('category')} (conf {label['data'].get('confidence')})")
The first label call for each feature is slow (~30s, ~10 Forge credits as the
labeler runs SAE on a 10-protein panel); subsequent calls hit cache.
Step 4-alt: permutation-test the SAE features (optional, more rigorous)
defpermutation_pvalues(cluster_positions, max_drop, n_perm=10000, rng=None):
"""Per-feature: is the cluster's mean drop > random equally-sized set?"""
rng = rng or np.random.default_rng(0)
n_positions = max_drop.shape[0]
cluster_size = len(cluster_positions)
observed = max_drop[cluster_positions].mean(axis=0)
null_geq = np.zeros(max_drop.shape[1], dtype=np.int32)
for _ inrange(n_perm):
idx = rng.choice(n_positions, size=cluster_size, replace=False)
null_geq += (max_drop[idx].mean(axis=0) >= observed).astype(np.int32)
return (null_geq + 1) / (n_perm + 1)
from statsmodels.stats.multitest import multipletests
p_raw = permutation_pvalues(np.array(cluster_cols), max_drop_per_pos)
_, p_adj, _, _ = multipletests(p_raw, method="fdr_bh")
significant_features = np.where(p_adj < 0.05)[0].tolist()
Use the mean of max_drop, not the max โ under this null a maximum-based
statistic returns almost no significant features. Single-position clusters
return no significant features (no statistical power) โ fall back to Step 4
descriptive ranking.
Step 5: Synthesize the mechanism call
For each cluster, combine the 3 (or 4 with SAE) evidence streams:
defcall_mechanism(structural, uniprot, sae_labels=None):
"""Return one of: catalytic | ligand-binding | interface |
structural-core | PTM | regulatory | mixed | unknown."""# Direct UniProt evidence winsif uniprot["active_site"]:
return"catalytic"if uniprot["binding_site"]:
return"ligand-binding"if uniprot["modified_residue"] andlen(uniprot["modified_residue"]) >= len(cluster) // 2:
return"PTM"# Structural evidenceif structural["ligand_pocket_count"] >= len(cluster) // 2:
return"ligand-binding"if structural["interface_count"] >= len(cluster) // 2:
return"interface"if structural["core_count"] >= len(cluster) // 2:
return"structural-core"# SAE evidence as tiebreakerif sae_labels:
from collections import Counter
cat_counts = Counter(l for l in sae_labels if l)
if cat_counts:
top_cat, n = cat_counts.most_common(1)[0]
if n >= 2: # at least 2 of top-5 SAE features agreereturn top_cat # e.g. "ligand-binding"return"unknown"
Step 6: Report
Cluster 1 โ positions [12, 13]
Structural: 0/2 interface, 2/2 ligand pocket (GTP), 0/2 core, all in P-loop helix
UniProt: Binding site (GTP) at residues 12, 13
Domain: small GTPase
SAE top 5: ligand-binding (ร2), secondary-structure (ร3)
โ MECHANISM: ligand-binding (GTP P-loop)
Cluster 2 โ positions [40, 41]
Structural: 2/2 interface (chain B = RAF1-RBD)
UniProt: no specific annotation
SAE top 5: structural-stability (ร3), domain (ร2)
โ MECHANISM: interface (KRAS-RAF1 binding)
Step 7: Visualize โ annotated DMS heatmap with hotspot callouts
The publication-style figure: DMS effect heatmap, sequence strip, structural
annotation track, with per-hotspot mechanism callouts above the heatmap.
Align everything to one position axis. Heatmap column p, sequence
letter p, every annotation bar covering residue p โ all share x = p.
Verify a landmark before drawing:
landmark_col = positions.index(12) # column for KRAS pos 12assert sequence[landmark_col] == "G", f"alignment broken at col {landmark_col}"
A 1-2 residue misalignment between heatmap and annotation track is a common,
visually subtle error. If you've cross-joined two coordinate systems and any
join was off-by-N, the whole figure is silently wrong. Verify here.
Heatmap + sequence + annotation track + callouts:
import matplotlib.pyplot as plt
import numpy as np
vlim = max(abs(np.nanmin(dms_matrix)), abs(np.nanmax(dms_matrix)))
fig, axes = plt.subplots(
nrows=4, ncols=1, figsize=(max(8, 0.15 * len(positions)), 6),
gridspec_kw={"height_ratios": [0.5, 4, 0.3, 0.5]}, sharex=True,
)
ax_callouts, ax_heat, ax_seq, ax_anno = axes
# Heatmap โ symmetric diverging (RdBu_r); center on 0 for ฮฮG-style data
im = ax_heat.imshow(
dms_matrix, aspect="auto", cmap="RdBu_r",
vmin=-vlim, vmax=vlim,
extent=(0, len(positions), 20, 0),
)
ax_heat.set_yticks(np.arange(20) + 0.5)
ax_heat.set_yticklabels(list(amino_acid_order))
ax_heat.set_ylabel("Substitution")
# Mark WT cells (box, no fill) โ distinguish "WT" from "not measured"for col, p inenumerate(positions):
wt_aa = sequence[col]
if wt_aa in amino_acid_order:
row = amino_acid_order.index(wt_aa)
ax_heat.add_patch(plt.Rectangle(
(col, row), 1, 1, fill=False, edgecolor='black', linewidth=0.5,
))
# Sequence strip โ one monospace letter per column
ax_seq.set_xlim(0, len(positions))
ax_seq.set_ylim(0, 1)
ax_seq.set_yticks([])
for col, letter inenumerate(sequence):
ax_seq.text(col + 0.5, 0.5, letter, ha="center", va="center",
family="monospace", fontsize=8)
# Annotation track โ region colors (top half) + core bar (bottom half)
anno_by_pos = {a["position"]: a for a in struct["data"]["annotations"]}
region_colors = {"interface": "#1f77b4", "ligand": "#ff7f0e",
"both": "#2ca02c", "other": "#cccccc"}
for col, p inenumerate(positions):
a = anno_by_pos.get(p, {})
ax_anno.add_patch(plt.Rectangle(
(col, 0.5), 1, 0.5, facecolor=region_colors.get(a.get("region", "other"), "#cccccc"),
))
if a.get("is_core"):
ax_anno.add_patch(plt.Rectangle((col, 0.0), 1, 0.5, facecolor="black"))
ax_anno.set_xlim(0, len(positions))
ax_anno.set_ylim(0, 1)
ax_anno.set_yticks([0.25, 0.75])
ax_anno.set_yticklabels(["core", "region"])
ax_anno.set_xlabel("Residue position")
# Callout row โ per-hotspot mechanism boxes linked to clusters by bracketsfor cluster, mechanism, top_features in hotspot_results:
cluster_cols = [positions.index(p) for p in cluster if p in positions]
ifnot cluster_cols:
continue
c_left, c_right = min(cluster_cols), max(cluster_cols)
center = (c_left + c_right) / 2
ax_heat.plot([c_left, c_right + 1], [0, 0], "k-", lw=2)
label_lines = [f"MECHANISM: {mechanism}"] + [f" {fl}"for fl in top_features[:3]]
ax_callouts.text(center, 0.5, "\n".join(label_lines),
ha="center", va="center", fontsize=7,
bbox=dict(facecolor="white", edgecolor="black"))
ax_callouts.plot([center, center], [0, -0.3], "k-", lw=0.5)
ax_callouts.set_xlim(0, len(positions))
ax_callouts.set_ylim(0, 1)
ax_callouts.axis("off")
fig.colorbar(im, ax=ax_heat, label="DMS effect (ฮฮG kcal/mol)")
plt.savefig("dms_hotspots_annotated.png", dpi=200, bbox_inches="tight")
Three cell-color rules to get right:
Real measurement โ diverging colour
WT cell โ boxed (the black outline above), value-cell colour = centre
Not measured โ distinct colour (e.g. light grey, not white โ white reads
as "neutral" against the diverging palette)
Long proteins: for >300 residues, split into multiple horizontal panels
(one panel per domain) rather than shrinking column width โ the per-residue
detail disappears below ~3 pixels per column.
Reproducing a published panel: verify its track alignment before
treating it as ground truth. Published DMS panels do carry registration
errors (the KRAS Fig 1i in the original paper is shifted +2 relative to its
own sequence โ see tooluniverse-protein-structural-annotation-pdb pitfalls).
Interpretation table โ what the mechanism call means downstream
Mechanism
Implication
catalytic
Direct enzyme function โ mutations abolish activity
ligand-binding
Substrate / cofactor / ion / nucleotide binding โ mutations alter substrate specificity or affinity
interface
Protein-protein interaction surface โ mutations may disrupt complex formation (consider PPI inhibitor design)
structural-core
Fold stability โ mutations destabilize protein (consider rescuing with chaperones; harder to drug)
PTM
Regulation site (phospho, acetyl, ubiquitin, glycosylation) โ mutations alter signaling rather than activity
regulatory
Allosteric site / autoinhibitory residue โ mutations bias conformational equilibrium
No mechanism could be assigned โ possibly novel function or wrong reference structure
Honest limitations
Wrong PDB โ wrong call. If your PDB doesn't include the relevant ligand
or partner, the structural evidence layer is blind to that mechanism. Pick
the structure that contains the right complex.
UniProt annotations are sparse for non-model proteins. Active sites are
well-curated for canonical enzymes; novel proteins may have no annotated
features and the skill falls back to structural + SAE only.
Single-position clusters limit statistical evidence. Descriptive
ranking still works but permutation p-values can't (n=1).
SAE feature labels are interpretive hints, not ground truth. Labels
come from how features activate across UniRef90, not per-protein expert
curation. Treat "category: ligand-binding" as a hypothesis weight, not a
proof.
Hotspots โ druggable sites. A catalytic residue is a critical residue
but not necessarily a good drug target (allosteric pockets often are
better). This skill explains why a residue is critical, not whether it's a
good target.
The mechanism call is a synthesis of evidence, not a measurement.
Don't quote the category as a fact โ quote the evidence and the call as a
reasoned conclusion.
