tooluniverse-residue-functional-mechanism-interpretation

残基功能机制解读

核心用户问题：“这些残基为何在功能上至关重要？”

回答这个问题需要不止一种来源：配体口袋中的残基与催化三联体、界面、疏水核心或 PTM 基序中的残基，意义各不相同。本技能综合多个 TU 工具提供的证据，为每个残基（或相邻残基簇）判定其机制。

残基可以来自任何来源——本技能不依赖残基的来源：

残基来源	常见调用方式
DMS 图谱热点	使用步骤 1（可选）按最大效应检测前 K 个热点，然后继续
ClinVar 复发变异	从 ClinVar 提取复发位点，直接作为 user_provided_positions 传入
文献热点区域	从论文图 1 中复制位点，直接传入
进化保守残基	按保守性评分筛选，传入前 N 个位点
可成药位点残基	从结合位点预测器获取，直接传入
临床医生的问题	“为什么 R175 位点的突变总是在肿瘤中出现？”——传入 [175]

---

何时使用此技能

• 你有一组残基（来源不限），需要解释它们在生物学上为何重要
• 你正在撰写论文的方法或讨论部分，需要逐个残基或残基簇提出机制性论断
• 你正在比较不同直系同源蛋白的残基集合，需要逐个残基的类别以便对齐

不适用于：

• 针对 DMS 验证预测器——使用 tooluniverse-variant-predictor-dms-validation
• 单个变异的 SAE 特征分解，用于查看哪些 ESMC 特征发生了变化——使用 tooluniverse-protein-sae-variant-interpretation
• 针对单个变异的 LoF 机制综合分析——使用 tooluniverse-protein-lof-mechanism

单变异技能侧重于一种突变的特征；本技能侧重于哪些残基重要，以及为何重要。

---

必需输入（选择一种入口路径）

路径 A — 直接提供残基（适用于 ClinVar、文献或 自定义位点）：

输入	说明
user_provided_positions: List[int]	从 1 开始计数的标准残基位点
蛋白质元数据	UniProt 登录号 + PDB ID + 链
（可选）DMS 矩阵	若提供，则用于为每个残基簇补充效应量背景
（可选）SAE 张量	若提供，则提供第 4 层 SAE 特征证据

路径 B — 从 DMS 图谱中检测热点（原始用例）：

输入	来源	说明
DMS 效应矩阵 (20, n_positions)	MaveDB_get_effect_matrix	未测量的位置为 NaN
disruptive_tail	DMS 检索元数据	"top" 或 "bottom"
蛋白质元数据	UniProt 登录号 + PDB ID + 链	用于多证据检索
（可选）SAE 证据	针对一个连续簇调用 ESM_get_region_sae_features（1 次 Forge 调用），或使用针对每个突变预先计算的完整 DMS SAE 张量（来自 ESM_get_sae_features）	SAE 证据层；参见步骤 4，了解应选择哪种路径

在路径 B 中，本技能会运行步骤 1 来检测热点；在路径 A 中，则会完全跳过步骤 1，直接进入步骤 2（收集证据）。

---

工作流

步骤 0（如果用户指定了具体位置，则为强制步骤）：前提核查

如果用户说“解释为什么残基／簇 X 是热点”，不要直接将此视为事实。请先核实它在本次 DMS 中是否确实是热点——用户可能会把其他情境中的生物学知识带入，而这些知识未必适用于特定实验测量的内容。

# 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")

判断规则：

• 如果指定位置按最大效应值排名位于前 25%，则前提得到确认，继续分析。
• 如果指定位置位于前 50% 但不在前 25%，则前提得到较弱支持；继续分析，但在报告中注明排名。
• 如果指定位置排名低于前 50%，则在继续机制分析之前，务必在回答开头透明地说明这一不符之处。

具体示例：用户问“为什么 KRAS G12/G13 在本次 DMS 中是折叠热点？”数据表明，按最大 ΔΔG 排名，G12 为 105/187（前 56%），G13 为 124/187（前 66%）——在本次 AbundancePCA 实验中，它们并非折叠热点（尽管它们确实是著名的致癌位置）。此技能的作用是先明确指出这一矛盾，然后继续分析这些残基的作用机制（它们的致癌效应是通过消除 GTPase 活性，而非破坏折叠实现的——这确实能有效回答用户真正关心的科学问题，只是与其字面提出的问题不同）。

步骤 1（仅适用于路径 B）：从 DMS 矩阵中检测热点

如果用户已提供位置（路径 A），则完全跳过此步骤。 在路径 A 中，残基列表就是输入；直接跳至本节底部的聚类步骤。

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 工作流约定（可用与不可用的内容）：

证据层	路径 A（用户提供的残基）	路径 B（DMS 热点）
结构信息（步骤 2）	✓ 始终可用	✓ 始终可用
UniProt 特征（步骤 3）	✓ 始终可用	✓ 始终可用
SAE 每特征标签（步骤 4，描述性）	✓ 如果提供 SAE 张量	✓ 如果提供 SAE 张量
SAE 置换检验（步骤 4-alt）	✗ — 需要从 DMS 得出的 max_drop 基线；对于没有 DMS 背景的位置没有意义	✓ 如果提供 SAE 张量
DMS 效应量背景	✓ 如果提供 DMS 矩阵（仅用于富集分析）	✓ 始终可用
机制综合分析（步骤 5）	✓ 始终可用	✓ 始终可用

步骤 2：收集每个簇的结构证据

对蛋白质结构进行一次注释，然后读取各簇位置对应的字段：

# 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_core")),
        "ss_elements": [by_pos.get(p, {}).get("ss_element") for p in cluster],
    }

步骤 3：收集每个簇的 UniProt 特征证据

# 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 in range(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 if any(f["type"] == "Active site" for f in features_by_pos.get(p, []))],
        "binding_site": [p for p in cluster if any(f["type"] == "Binding site" for f in features_by_pos.get(p, []))],
        "modified_residue": [p for p in cluster if any(f["type"] == "Modified residue" for f in features_by_pos.get(p, []))],
        "disulfide_bond": [p for p in cluster if any(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"],
    }

步骤 4（可选）：对每个簇进行 SAE 特征排名

根据是否已有完整的 DMS SAE 张量，采用以下两种路径：

路径 A — 没有预计算张量（大多数情况）：直接使用区域工具。 如果簇是连续范围（或可以扩展为连续范围），ESM_get_region_sae_features 会在一次 Forge 调用中汇总该范围内的 SAE 特征：

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"]]

这是推荐的默认方案——与 DMS 张量路径相比，只需 1 次 Forge 调用，而非 20 × 簇大小。对于非连续簇，请针对每个连续子范围分别运行一次，并合并排名前 K 的特征。

路径 B — 已有来自 DMS 扫描的预计算 SAE 张量（例如，variant-predictor-dms-validation 管道在磁盘上留下了一个张量）：直接计算下降值，无需再次调用 Forge。

def cluster_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)

使用 SAE 特征标注器为每个排名靠前的特征添加标签：

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')})")

每个特征的首次标注调用速度较慢（约 30 秒，约消耗 10 个 Forge 点数，因为标注器会在一个包含 10 种蛋白质的面板上运行 SAE）；后续调用会命中缓存。

步骤 4-alt：对 SAE 特征进行置换检验（可选，更严谨）

def permutation_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 _ in range(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()

使用 max_drop 的均值，而非最大值——在此零假设下，基于最大值的统计量几乎不会产生显著特征。单位置簇不会返回显著特征（没有统计效力）——改用步骤 4 的描述性排名。

步骤 5：综合判断作用机制

对每个簇，综合 3 条证据流（如果包含 SAE，则为 4 条）：

def call_mechanism(structural, uniprot, sae_labels=None):
    """Return one of: catalytic | ligand-binding | interface |
    structural-core | PTM | regulatory | mixed | unknown."""

    # Direct UniProt evidence wins
    if uniprot["active_site"]:
        return "catalytic"
    if uniprot["binding_site"]:
        return "ligand-binding"
    if uniprot["modified_residue"] and len(uniprot["modified_residue"]) >= len(cluster) // 2:
        return "PTM"

    # Structural evidence
    if 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 tiebreaker
    if 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 agree
                return top_cat  # e.g. "ligand-binding"

    return "unknown"

步骤 6：报告

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)

第 7 步：可视化 — 带热点标注的 DMS 热图

这是一张适合发表的图：DMS 效应热图、序列条、结构注释轨道，以及热图上方针对各个热点的机制标注。

将所有内容对齐到同一个位置坐标轴。 热图列 p、序列字母 p，以及覆盖残基 p 的每个注释条，都共享 x = p。绘制前先验证一个标志位点：

landmark_col = positions.index(12)  # column for KRAS pos 12
assert sequence[landmark_col] == "G", f"alignment broken at col {landmark_col}"

热图和注释轨道之间错位 1–2 个残基，是一种常见且视觉上不易察觉的错误。如果你对两个坐标系做了交叉连接，而且任何一次连接都偏移了 N 位，那么整张图就会悄无声息地出错。请在此验证。

热图 + 序列 + 注释轨道 + 标注：

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 in enumerate(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 in enumerate(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 in enumerate(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 brackets
for cluster, mechanism, top_features in hotspot_results:
    cluster_cols = [positions.index(p) for p in cluster if p in positions]
    if not 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")

需要正确处理的三条单元格着色规则：

• 实际测量值 → 发散色
• WT 单元格 → 加框（如上方的黑色轮廓），数值单元格的颜色 = 中心色
• 未测量 → 使用独特颜色（例如浅灰色，而非白色——在发散色板中，白色会被读作“中性”）

长蛋白质：对于超过 300 个残基的蛋白质，应拆分成多个水平面板（每个结构域一个面板），而不是缩小列宽——每残基的细节在每列小于约 3 像素时就会消失。

复现已发表的面板：将其轨道对齐视为准确依据之前，先验证对齐情况。已发表的 DMS 面板确实存在配准错误（原论文中的 KRAS Fig 1i 相对于其自身序列偏移了 +2——参见 tooluniverse-protein-structural-annotation-pdb 的注意事项）。

---

解读表 — 机制判断对后续工作的含义

机制	含义
催化	直接影响酶功能——突变会消除活性
配体结合	底物 / 辅因子 / 离子 / 核苷酸结合——突变会改变底物特异性或亲和力
界面	蛋白质-蛋白质相互作用表面——突变可能破坏复合物形成（可考虑设计 PPI 抑制剂）
结构核心	折叠稳定性——突变会使蛋白质不稳定（可考虑用伴侣蛋白挽救；更难开发药物）
PTM	调控位点（磷酸化、乙酰化、泛素化、糖基化）——突变改变的是信号传导，而非活性
调控	变构位点 / 自抑制残基——突变会改变构象平衡
混合	多种证据类型相互矛盾——需要逐个案例分析
未知	无法判定机制——可能是新功能，也可能是参考结构不正确

---

诚实说明局限

1. PDB 不正确 → 判断不正确。如果你的 PDB 不包含相关配体或伙伴蛋白，结构证据层就无法识别该机制。请选择包含正确复合物的结构。
2. 对于非模式蛋白，UniProt 注释较少。典型酶的活性位点经过充分整理；新型蛋白可能没有已注释特征，此时该技能只能依赖结构证据 + SAE。
3. 单位置簇会限制统计证据。描述性排序仍然有效，但无法计算置换 p 值（n=1）。
4. SAE 特征标签是用于解读的提示，不是既定事实。这些标签依据特征在 UniRef90 中的激活情况生成，而非由专家针对每种蛋白质逐一整理。应将“category: ligand-binding”视为假设的权重，而非证明。
5. 热点 ≠ 可成药位点。催化残基是关键残基，但不一定是理想的药物靶点（变构口袋往往更好）。此技能解释残基为何关键，而不判断它是否是理想靶点。
6. 机制判断是对证据的综合，而非测量结果。不要将类别作为事实引用——应引用证据，并将判断表述为有理有据的结论。

---

交叉引用

工具 / 技能	作用
MaveDB_get_effect_matrix	DMS 矩阵输入
tooluniverse-protein-structural-annotation-pdb（或直接使用 Structure_annotate_per_residue）	结构证据
UniProt_get_function_by_accession	UniProt 特征（活性位点、结合位点、PTM、二硫键）
ESM_get_region_sae_features	步骤 4 路径 A — 在 1 次 Forge 调用中聚合连续簇上的 SAE 特征（首选）
ESM_get_sae_features	步骤 4 路径 B — 仅当你已有预计算的完整 DMS SAE 张量时使用
ESM_describe_sae_feature	在步骤 4 中标注 SAE 特征
tooluniverse-variant-predictor-dms-validation	相关技能：在信任预测结果之前验证预测器
（热图可视化现为此技能的步骤 7）	带有逐热点标注的 DMS 面板
alphafold_get_prediction	没有实验 PDB 时提供 pLDDT 背景信息