tooluniverse-residue-functional-mechanism-interpretation

残基の機能的なメカニズムの解釈

ユーザーが最も知りたいこと：「これらの残基が機能上重要なのはなぜですか？」

答えるには、1つの情報源だけでは不十分です。リガンドポケット内の残基と、触媒三残基、相互作用面、疎水性コア、PTM モチーフ内の残基とでは、意味が異なります。このスキルでは、複数の TU ツールから得られる証拠を統合し、各残基（または隣接する残基のクラスター）のメカニズムを特定します。

残基の出所は問いません。このスキルは、その出所に依存しません。

残基の出所	典型的な呼び出しパターン
DMS マップのホットスポット	Step 1（任意）を使って最大効果の上位 K 件を検出し、その後続行します
ClinVar の再発性バリアント	ClinVar から再発位置を取得し、user_provided_positions として直接渡します
文献に記載されたホット領域	論文の Fig 1 から位置を貼り付けて、直接渡します
進化的に保存された残基	保存スコアで絞り込み、上位 N 件の位置を渡します
創薬可能部位の残基	結合部位予測ツールから取得し、直接渡します
臨床医からの質問	「腫瘍で R175 の変異が繰り返し見られるのはなぜですか？」— [175] を渡します

---

このスキルを使用する場面

• どこから得たものでも、残基のリストがあり、生物学的になぜ重要なのか説明する必要がある場合
• 論文の方法または考察のセクションを執筆しており、残基またはクラスターごとのメカニズムに関する主張が必要な場合
• オーソログ間で残基セットを比較しており、対応付けに使う残基ごとのカテゴリが必要な場合

対象外：

• DMS を使って予測器を検証する場合 — tooluniverse-variant-predictor-dms-validation を使用します
• ESMC のどの特徴が変化したかを確認したい単一バリアントのSAE 特徴分解 — tooluniverse-protein-sae-variant-interpretation を使用します
• 単一バリアントについて、単独でLoF メカニズムを統合する場合 — tooluniverse-protein-lof-mechanism を使用します

単一バリアント向けのスキルは1つの変異の特徴に焦点を当てます。このスキルは、どの残基が、なぜ重要なのかに焦点を当てます。

---

必須入力（どちらかの入力経路を選択）

経路 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 の証拠	連続した1つのクラスターに対する ESM_get_region_sae_features（Forge 呼び出し1回）、または変異体ごとの ESM_get_sae_features による事前計算済みの DMS 全体の SAE テンソル	SAE の証拠レイヤー。どちらの経路を使うかは Step 4 を参照してください

経路 B では、このスキルは Step 1 を実行してホットスポットを検出します。経路 A では Step 1 を完全に省略し、Step 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％未満の場合 → この不一致を明確に報告してください。 メカニズム分析に進む前に、回答の冒頭で説明してください。

具体例：ユーザーが「この DMS では KRAS G12/G13 がなぜフォールディングのホットスポットなのですか？」と尋ねたとします。データでは、最大 ΔΔ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 テンソルがすでにあるかどうかに応じて、2 つの方法があります。

パス A — 事前計算済みのテンソルがない場合（ほとんどのケース）：領域ツールを直接使用します。 クラスターが連続した範囲である（または範囲をパディングできる）場合、ESM_get_region_sae_features を使うと、1 回の 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"]]

これが標準の方法です。Forge の呼び出しは 1 回で済みますが、DMS テンソルを使う方法では 20 × クラスターのサイズ分必要です。クラスターが連続していない場合は、連続した部分範囲ごとに 1 回実行し、上位 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')})")

ラベラーは 10 個のタンパク質のパネルで SAE を実行するため、各特徴についての最初のラベル付け呼び出しには時間がかかります（約 30 秒、Forge クレジット約 10）。その後の呼び出しはキャッシュが使われます。

ステップ 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 効果のヒートマップ、配列の帯、構造注釈トラックを配置し、ヒートマップの上にホットスポットごとのメカニズムを吹き出しで示します。

すべてを1つの位置軸に揃えてください。 ヒートマップの列 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残基ずれるのは、よくある目立ちにくいミスです。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")

セルの色に関する3つのルール：

• 実測値 → 発散色
• WT セル → 枠で囲む（上の黒い輪郭）。値セルの色は中心値にします
• 未測定 → 他と区別できる色（例：白ではなく薄いグレー。発散色パレットの中では、白は「中立」と受け取られるためです）

長いタンパク質：残基数が300を超える場合、列幅を縮めるのではなく、複数の横長パネル（ドメインごとに1パネル）に分割してください。1列あたり約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 — 連続クラスター全体の SAE 特徴を 1 回の Forge 呼び出しで集約します（推奨）
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 の参考情報