tooluniverse-variant-predictor-dms-validation

DMS に対するバリアント効果予測ツールのベンチマーク

ユーザーが抱く中心的な疑問：「バリアント効果予測ツールを使っていますが、このタンパク質の実験的な DMS 測定値と実際に相関していますか？」

予測ツールは、単一のミスセンスバリアントに数値スコアを割り当てるものであれば、どのようなものでもかまいません。

• 変異部位における ESM-C 6B Sparse Autoencoder（SAE）の特徴量の低下
• AlphaMissense の病原性スコア
• ESM ロジットベースのバリアントスコアリング（ESM_score_sequence）
• EVE / EVE++ スコア
• 保存性スコア（ConSurf、Rate4Site）
• DynaMut2 の ΔΔG 予測
• 独自の社内モデル

このスキルでは、同じ統計フレームワークを使って、これらすべてを DMS データセットに対して検証します。このコレクションの関連スキルが SAE をテーマにしているため、実例には SAE を取り上げていますが、手順は予測ツールに依存しません。

---

このスキルを使う場面

• バリアント効果予測ツールを選び、関心のあるタンパク質に対して信頼できるか確認したい場合
• 同じ DMS データセットで 2 つの予測ツールを比較する場合（このスキルを 2 回実行し、得られた K ごとの Mann-Whitney U 検定の p 値と効果量を比較してください）
• レビュアーから頑健性の根拠を求められている場合 — パラメーターのスイープ（K、中立帯、破壊的変異の分位点）がその根拠になります
• 新しいバリアント効果予測手法を公開するにあたり、ベンチマーク図が必要な場合

対象外：

• 位置ごと／特徴量ごとの解釈 — tooluniverse-residue-functional-mechanism-interpretation を使用してください
• 単一バリアントの解釈（バリアントが 1 つだけで、DMS がない場合）— tooluniverse-protein-sae-variant-interpretation または tooluniverse-protein-lof-mechanism を使用してください
• 予測ツールをゼロから構築すること — 対象外です

---

必要な入力

入力	形式	例
DMS 効果行列	（20 個のアミノ酸 × n_positions）np.array、未測定箇所は NaN	MaveDB_get_effect_matrix の出力
破壊的変異の裾の規約	"top"（ΔΔG が正 = 不安定化）または "bottom"（フィットネスが低い = LoF）	DMS 取得ステップのメタデータ
バリアントごとの予測スコア	（20 × n_positions）DMS と同じレイアウトの np.array	選択した予測ツールで算出 — ステップ 2 を参照
集約 K	int。予測ツールがバリアントごとに複数のサブスコアを出力する場合、上位 K 個の平均を取ります	SAE のような複数特徴量の予測ツールにのみ関係します。それ以外では無視してください

---

ワークフロー

ステップ 1：DMS を効果行列として取得する

1 回の呼び出しで、解析可能な (20 × n_positions) 行列が返されます — HGVS の解析、単一ミスセンス変異のフィルタリング、スコアフィールドの検出、（任意の）UniProt 番号付けの検証はツール内で行われます：

r = MaveDB_get_effect_matrix(
    urn="urn:mavedb:00000115-a-7",
    uniprot_id="P01116",     # optional but recommended — enables numbering check
)
matrix = np.array(r["data"]["matrix"], dtype=np.float32)  # (20, n_positions)
positions = r["data"]["positions"]
amino_acid_order = r["data"]["amino_acid_order"]          # always 'ACDEFGHIKLMNPQRSTVWY'
# Audit fields — surface in your report:
#   r["data"]["n_parsed_single_missense"], n_dropped, score_field_used,
#   numbering_offset, numbering_check

numbering_offset != 0 の場合、MaveDB の位置番号は UniProt の標準配列と異なります — 他のソース（PDB / SAE 特徴量 / AlphaMissense）と結合する前に、オフセットを適用してください。

ステップ 2：バリアントごとの予測スコアを算出する

次の予測元から 1 つ選んでください。選択によって変わるのはステップ 2 のみです。ステップ 3～5 は共通です。

Predictor option A — ESM-C 6B SAE（実例）

各バリアントについて、残基ウィンドウ全体のSAEアクティベーションを合計し、 drop = max(0, WT − mut) を計算して、dropの上位K個の平均により バリアントごとに1つのスコアに集約します。

import numpy as np

def sae_drop_per_variant(wt_pooled, variant_pooled, K=3):
    """SAE drop for one variant = mean of the K largest feature drops."""
    drops = np.maximum(0.0, wt_pooled - variant_pooled)  # (n_features,)
    sorted_desc = -np.sort(-drops)
    return float(sorted_desc[:K].mean())

飽和スイープのスコアリング方法は2つ — バッチサイズに応じて選択します。

条件	使用するもの	Forgeコスト（19種類の代替 × N位置）
バリアントが100個以下の飽和スイープ、またはバリアントごとの上位K個の差分のみが必要な場合	ESM_score_variant_sae_batch(sequence, variants=[...], top_k_features=10)	1 + 19N（1参照 + 変異ごとに1）
タンパク質全体の特徴量テンソル（後続のPCA／クラスタリングなどに使用）	ESM_get_sae_features(sequence=mutant) をループで実行し、(sequence, position) でキャッシュ	2 × 19N（キャッシュ済みの再実行は無料）

バッチツールをデフォルトとして使用するのが適切です。バリアントごとに機能を調べるパターンと比べてForgeコストを半分に抑えられます。また、1回の呼び出しあたり100バリアントという上限は、1位置での飽和変異導入（19種類）や、短い位置スイープ（例：位置10～15では90バリアント）をカバーします。より長いスイープでは、複数のバッチ呼び出しに分割してください。

一部の予測解析で必要となる、残基ごと × 特徴量ごとの完全なテンソルには、ループパターンを使用してください。ライブラリ全体の規模でも十分にテストされています。tests/integration/test_dms_pipeline_e2e_kras.py では、KRASの位置10～25について約300個の変異体を実行します。

Predictor option B — AlphaMissense（hegelabプロキシ：置換ごとの数値ではなくカテゴリ値）

TU AlphaMissenseツールは、置換ごと（位置、alt_aa）の数値スコアではなく、残基レベルのカテゴリ割り当てを返す公開API（hegelab.org）へのプロキシです。利用できる呼び出しは3つあります。規模に合ったものを選択してください。

ツール	シグネチャ	返される値
AlphaMissense_get_variant_score	uniprot_id、variant（例："p.G12V"）	alt_aaのビン（benign／ambiguous／pathogenic）を含む残基レベルのデータ。また、mean／mean_allも返します
AlphaMissense_get_residue_scores	uniprot_id、position	同じ残基レベルのデータ（位置ごとの検索で、バリアント呼び出しを19回行うより低コストです）
AlphaMissense_get_protein_scores	uniprot_id	1回の呼び出しでタンパク質全体のデータを取得します（タンパク質全体のDMS解析では最も低コストです）

3つの呼び出しすべてで、レスポンスの形式は同じです。KRASの12番目の位置の例：

r = AlphaMissense_get_residue_scores(uniprot_id="P01116", position=12)
# r["data"]["scores"]:
#   {"uid":"P01116","aa":"G","resi":12,
#    "benign":"", "ambiguous":"",
#    "pathogenic":"6:A,C,D,R,S,V",                  # ← SNV-reachable subs only
#    "pathogenic_all":"19:A,C,D,E,F,...,Y",         # ← all 19 substitutions
#    "mean":0.9885,           # mean over SNV-reachable subs
#    "mean_all":0.9950}       # mean over all 19 substitutions
# r["data"]["thresholds"]:
#   {"pathogenic":"> 0.564","ambiguous":"0.34 - 0.564","benign":"< 0.34"}

(20, n_positions) 予測行列に値を入力するには、ビン文字列を解析し、各alt_aaを数値スコアにマッピングします（ツールでは置換ごとの真の数値スコアが公開されていないため、標準的な補完方法としてビンの中点を使用します）。

BIN_MIDPOINTS = {"benign": 0.17, "ambiguous": 0.452, "pathogenic": 0.782}

def parse_bin_list(bin_str):
    """'6:A,C,D,R,S,V' → ['A','C','D','R','S','V']; '' → []."""
    if not bin_str:
        return []
    _, aas = bin_str.split(":", 1)
    return aas.split(",")

def am_per_variant_matrix(uniprot_id, positions, aa_index):
    AAS = list(aa_index)
    M = np.full((20, len(positions)), np.nan, dtype=np.float32)
    # One whole-protein call is much cheaper than n_positions calls
    p = AlphaMissense_get_protein_scores(uniprot_id=uniprot_id)
    per_res = {row["resi"]: row for row in p["data"]["scores"]}
    for pos_idx, pos in enumerate(positions):
        row = per_res.get(pos, {})
        for cat in ("benign", "ambiguous", "pathogenic"):
            for alt in parse_bin_list(row.get(f"{cat}_all", "")):
                if alt in aa_index:
                    M[aa_index[alt], pos_idx] = BIN_MIDPOINTS[cat]
    return M

より高解像度の代替手段 — DeepMindの一括CSV（ビンの中点ではなく、置換ごとの真の数値がどうしても必要な場合に限ります）。

# The official DeepMind release contains per-variant continuous scores
# (not just bin assignments). Not currently wrapped by a TU tool — fetch directly:
#   https://alphafold.ebi.ac.uk/files/AF-<uniprot_id>-F1-aa-substitutions.csv
# or stream-filter https://storage.googleapis.com/dm_alphamissense/AlphaMissense_aa_substitutions.tsv.gz

順位相関の解析が必要な場合は、一括CSVを使用してください（Spearman相関には、3つのビンの中点ではなく連続値が適しています）。「このバリアントは病原性ビンに含まれるか」という二値の情報だけが必要な場合は、TUプロキシを使用してください。

どちらの方法も無料で、APIキーは不要です。

Predictor option C — ESMロジットベースのスコア

ESM_score_sequence(
    sequence=mutant_sequence,
    model="esmc-600m-2024-12",
)
# returns per-residue logits; compute mutant-vs-WT log-odds at the mutation site

Predictor option C2 — ESM-2 masked-marginal（キー不要、APIキー不要）

ESM_API_KEY がない場合（オプションA／Cのツールにはキーが必要です）、ESM2_score_missense_variant を使用してください。HuggingFaceの無料プロバイダーhf-inference経由でESM-2を実行し、1つのミスセンスバリアントについてmasked-marginal対数尤度比logP(mut) − logP(wt)（Meier et al. 2021）を返します。

# one call per variant; negative LLR = mutant disfavored (candidate deleterious)
res = ESM2_score_missense_variant(
    sequence=wild_type_sequence,   # 1-letter AA string
    position=position,             # 1-based
    mutant=alt_aa,                 # e.g. "V"
)
score = res["data"]["log_likelihood_ratio"]   # use directly as the predictor score

バリアントごとにHTTP呼び出しを1回行う必要があり（一括エンドポイントはありません）、飽和スイープではキーを使うESM-Cのバッチツールよりも時間がかかります。キーがある場合はオプションA／Cを優先し、セットアップ不要の代替手段としてこちらを利用してください。約1022残基を超える配列は、バリアントを中心としたウィンドウに自動的に分割されます（metadata.windowedを参照）。LLRは順位付け用のスコアであり、校正済みの確率ではありません。閾値処理はステップ3～5で行うため、生のLLRをそのまま(20, n_positions)行列に入力してください。

Predictor option D — 任意の外部スコア

独自のスコアをご用意ください。DMS行列に位置を揃えた(20, n_positions)のnp.ndarrayを作成してください。

ステップ 3：DMS を中立群と破壊群に層別化する

def categorize(dms_matrix, disruptive_tail, neutral_abs=0.1, disruptive_quantile=0.05):
    """Split variants into neutral and disruptive masks.

    disruptive_tail: 'top' (positive = destabilizing, e.g. folding ΔΔG)
                     'bottom' (low = LoF, e.g. fitness)
    """
    flat = dms_matrix[~np.isnan(dms_matrix)]
    if disruptive_tail == "top":
        cut = np.quantile(flat, 1 - disruptive_quantile)
        disruptive = dms_matrix >= cut
    elif disruptive_tail == "bottom":
        cut = np.quantile(flat, disruptive_quantile)
        disruptive = dms_matrix <= cut
    else:
        raise ValueError("disruptive_tail must be 'top' or 'bottom'")
    neutral = np.abs(dms_matrix) <= neutral_abs
    disruptive = disruptive & ~np.isnan(dms_matrix)
    neutral = neutral & ~np.isnan(dms_matrix) & ~disruptive
    return neutral, disruptive

中立帯は狭く保ってください（|effect| ≤ 0.1）。中立帯が広すぎると、弱い破壊的変異が「中立」群に混入し、対比が弱まります。 符号は重要です — DMS 検索スキルのメタデータから disruptive_tail を取得してください。符号が逆だと、すべての結論が気づかないうちに反転します。

ステップ 3.5（必須）：MWU 実行前の健全性チェック — NaN によるサイレントな失敗を検出する

統計検定を実行する前に、予測スコアが実際に入力されていることを確認してください。サイレントな NaN 行列は、このワークフローで最もよく発生する失敗モードです。途中でエラーになった SAE のバッチ計算、変異をスキップした AlphaMissense の取得、タイムアウトした ESM forge の呼び出しなどが原因となり、意味のない統計値を報告する「成功」した実行につながります。

neutral, disruptive = categorize(dms_matrix, disruptive_tail)
s_neutral_all = predictor_scores[neutral]
s_disruptive_all = predictor_scores[disruptive]
s_neutral_finite = s_neutral_all[~np.isnan(s_neutral_all)]
s_disruptive_finite = s_disruptive_all[~np.isnan(s_disruptive_all)]

# Hard gate: NaN coverage check
coverage_n = len(s_neutral_finite) / max(len(s_neutral_all), 1)
coverage_d = len(s_disruptive_finite) / max(len(s_disruptive_all), 1)
if len(s_neutral_finite) < 5 or len(s_disruptive_finite) < 5:
    raise ValueError(
        f"Insufficient predictor scores: only {len(s_neutral_finite)} neutral "
        f"and {len(s_disruptive_finite)} disruptive non-NaN values. "
        f"Coverage = {coverage_n:.0%} / {coverage_d:.0%}. "
        f"Predictor matrix may be empty / failed to populate. "
        f"INVESTIGATE THE PREDICTOR COMPUTATION STEP before continuing."
    )
if coverage_n < 0.5 or coverage_d < 0.5:
    print(f"WARNING: predictor coverage only {coverage_n:.0%} (neutral) / "
          f"{coverage_d:.0%} (disruptive). Results may be biased toward the "
          f"non-NaN subset.")

ValueError が発生した場合、「予測器 X がデフォルトで勝った」とごまかさないでください。予測器の計算手順をデバッグ（再実行、API キーの確認、ログファイルの確認）し、何が失敗したのかを報告するのが適切な対応です。

符号規約の再確認（データセットが新しい場合も必須です）：内部の基準点を使って、データに対する disruptive_tail を確認してください。メタデータのフィールドは、サブセットでは誤解を招くことがあります（例：TP53 DBD のウィンドウでは、TP53 全体の存在量アッセイと逆方向になる場合があります）。最も簡単な確認方法は、固定された方向に一致することがわかっている予測器（AlphaMissense の病原性スコアは常に正の値ほど損傷が大きい）と、DMS 効果の Spearman 相関を計算することです。

from scipy.stats import spearmanr
flat_dms = dms_matrix.ravel()
flat_pred = predictor_scores.ravel()
mask = ~(np.isnan(flat_dms) | np.isnan(flat_pred))
rho, p_corr = spearmanr(flat_dms[mask], flat_pred[mask])
expected_sign = "+" if disruptive_tail == "top" else "-"
got_sign = "+" if rho > 0 else "-"
if got_sign != expected_sign and abs(rho) > 0.1:
    print(f"WARNING: Spearman rho = {rho:.3f} (sign={got_sign}) but "
          f"disruptive_tail='{disruptive_tail}' expects sign={expected_sign}. "
          f"Sign convention may be inverted for THIS subset. Verify before "
          f"interpreting MWU.")

ステップ 4：片側 Mann-Whitney U 検定

from scipy.stats import mannwhitneyu

u, p = mannwhitneyu(s_disruptive_finite, s_neutral_finite, alternative="greater")
print(f"disruptive median = {np.median(s_disruptive_finite):.4f}, "
      f"neutral median = {np.median(s_neutral_finite):.4f}, p = {p:.3g}")

top-K パラメーターを持つ SAE スタイルのマルチ特徴予測器では、K ∈ {1, 3, 10} について実行し、3 つすべてを報告してください。最適な K は通常、予測器ごとに異なります。調整した結果ではなく、比較を透明にすることが重要です。

MWU は識別検定ですが、有効な分析方法はこれだけではありません。 次の方法の一方または両方も報告することを検討してください。

• 予測器と DMS 効果の Spearman 順位相関 — 二値の中立／破壊の分類では失われる段階的なキャリブレーションを捉えます。「この予測器はキャリブレーションされていますか？」という問いに最適で、「識別できますか？」だけを問う場合とは異なります。
• 臨床的に意味のある破壊の閾値（例：ΔΔG > 1 kcal/mol）における AUROC — 実務者に馴染みのある 0～1 の値で、文献にある予測器と直接比較できます。

このスキルのきっかけとなった評価（KRAS フォールディングの AlphaMissense ベンチマーク）では、MWU（p=0.20、「信頼性なし」）と Spearman（ρ=0.23、p<1e-30、「弱いキャリブレーション」）から、定性的には類似した結果が得られました。ただし、Spearman の解釈のほうが豊富な情報を含んでいました。ユーザーが単に「X は信頼できますか？」と尋ねた場合、データの形状に明らかに適さない方法がない限り、両方を提示してください。

ステップ 5：頑健性のスイープ

sweep = []
for neutral_abs in (0.05, 0.1, 0.2):
    for q in (0.05, 0.1):
        neut, disr = categorize(dms_matrix, disruptive_tail,
                                neutral_abs=neutral_abs,
                                disruptive_quantile=q)
        s_n = predictor_scores[neut][~np.isnan(predictor_scores[neut])]
        s_d = predictor_scores[disr][~np.isnan(predictor_scores[disr])]
        if len(s_n) < 5 or len(s_d) < 5:
            continue
        _u, p = mannwhitneyu(s_d, s_n, alternative="greater")
        sweep.append({"neutral_abs": neutral_abs, "disruptive_q": q, "p": p,
                      "n_n": len(s_n), "n_d": len(s_d)})

（neutral_abs、q）の一点でしか基準を満たさない予測器は疑わしいと考えられます。グリッド全体で基準を満たす予測器は頑健です。

ステップ 6：可視化して解釈する

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(4, 4))
ax.boxplot(
    [s_neutral, s_disruptive],
    labels=[f"neutral (n={len(s_neutral)})", f"disruptive (n={len(s_disruptive)})"]
)
ax.set_ylabel("Predictor score")
ax.set_title(f"p = {p:.3g}")
plt.tight_layout()
plt.savefig("predictor_vs_dms.png", dpi=150)

---

解釈

結果	意味
p < 0.01 かつスイープ全体で頑健	このタンパク質では、予測器が破壊的変異と中立変異を確実に識別します。優先順位付けには安全に使用できます（分類には使用できません。制限事項を参照）
p < 0.05 だがスイープ中に符号が反転	境界的です。効果は存在しますが、層別化に左右されます。より多くのデータ、またはより狭い中立帯が必要です
p > 0.05	この DMS アッセイでは予測器に情報量がありません。考えられる原因：disruptive_tail の誤り、タンパク質ファミリーに対する予測器のキャリブレーション不良、予測器の学習対象ではない特性を DMS が測定している
K=1 では有意だが、K=10 では有意でない（マルチ特徴予測器のみ）	破壊効果は少数の特徴に集中しています（K=1 が最適）。複数の特徴に分散している場合は K=10 が最適です

---

2つの予測器を比較する

同じ DMS データセット（例：SAE のドロップと AlphaMissense）に対して、このスキルを2回実行し、 同じ層別化（disruptive_tail、neutral_abs、disruptive_quantile）を維持して、 以下を比較します。

指標	より優れた予測器の特徴
より低い p 値	識別への確信度が高い
より大きい中央値の差（disruptive − neutral）	効果量が大きい
カバレッジ内の NaN が少ない	より多くのバリアントを予測する
スイープ全体で頑健	条件が異なっても信頼性が高い

---

正直に把握しておくべき限界

1. タンパク質ごとの検定のみです。「予測器 X は優れている」という一般化には、 異なるファミリーの複数のタンパク質で実行する必要があります。
2. アッセイごとの検定のみです。安定性向けに調整された予測器は、 結合適応度アッセイでは機能しない可能性があります。同じタンパク質でも、異なる DMS では結果が異なります。
3. シグナルは大まかなものです。「予測器が変異による破壊性に反応する」ことを示します。 「予測器が正しい残基を特定する」ことを示すものではありません。そのためには、 tooluniverse-residue-functional-mechanism-interpretation を実行してください。
4. MWU は独立性を前提としています。同じ位置にある DMS バリアントは、共有する構造的背景によりある程度相関するため、p 値はやや楽観的です。
5. テールの定義は恣意的です。破壊性のカットオフ 5% と中立バンド 0.1 は妥当なデフォルト値ですが、 適切な閾値はアッセイのシグナル対ノイズ比によって異なります。頑健性をもたらすのは、 単一の選択ではなくスイープです。

---

相互参照

ステップ	ツール／スキル
DMS の取得	MaveDB_get_effect_matrix
バリアントごとの SAE スコアリング（100 バリアント以下）	ESM_score_variant_sae_batch（推奨 — N+1 回の呼び出し）
バリアントごとの SAE スコアリング（テンソル全体／無制限）	ESM_get_sae_features（ループ＋キャッシュ）、ESM_score_variant_sae_disruption（単一バリアント）
バリアントごとの AlphaMissense（単一検索）	AlphaMissense_get_variant_score(uniprot_id, variant)
位置ごとの AlphaMissense（飽和変異導入）	AlphaMissense_get_residue_scores(uniprot_id, position)
タンパク質全体の AlphaMissense（DMS に最も低コスト）	AlphaMissense_get_protein_scores(uniprot_id)
バリアントごとの ESM ロジット	ESM_score_sequence
構造上の事前情報（予測器分析用）	Structure_annotate_per_residue
次のステップ：ホットスポットごとのメカニズム	tooluniverse-residue-functional-mechanism-interpretation
最終的な可視化	tooluniverse-residue-functional-mechanism-interpretation のステップ 7（注釈付きヒートマップ＋吹き出し）