tooluniverse-variant-predictor-dms-validation

针对 DMS 的变异效应预测器基准测试

核心用户问题：“我有一个变异效应预测器——它在这个蛋白质上的结果，真的与实验 DMS 测量值相关吗？”

预测器可以是任何一种能够为单个错义变异分配数值分数的工具：

• ESM-C 6B Sparse Autoencoder (SAE) 在突变位点的特征下降值
• AlphaMissense 致病性分数
• 基于 ESM logits 的变异评分（ESM_score_sequence）
• EVE / EVE++ 分数
• 保守性分数（ConSurf、Rate4Site）
• DynaMut2 ΔΔG 预测值
• 自定义的内部模型

此技能使用相同的统计框架，通过 DMS 数据集验证所有这些预测器。由于本系列中的相关技能都以 SAE 为主题，因此这里以 SAE 为示例演示；但此流程不依赖特定预测器。

---

何时使用此技能

• 你选择了一个变异效应预测器，想知道是否值得信任它在目标蛋白质上的预测结果
• 你要在同一个 DMS 数据集上比较两个预测器（运行此技能两次，然后比较每个 K 值对应的 Mann-Whitney U p 值和效应大小）
• 审稿人要求提供稳健性证据——参数扫描（K、中性区间、破坏性分位数）正好可以满足这一要求
• 你要发表一种新的变异效应方法，并需要一张基准测试图

不适用于：

• 按位点 / 按特征进行解读——请使用 tooluniverse-residue-functional-mechanism-interpretation
• 解读单个变异（你只有一个变异，没有 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 检索步骤的元数据
每个变异的预测器分数	与 DMS 布局匹配的 (20 × n_positions) np.array	使用所选预测器计算——见步骤 2
聚合 K	int，当预测器为每个变异输出多个子分数时，取前 K 个的均值	仅适用于 SAE 等多特征预测器；其他情况忽略

---

工作流

步骤 1：将 DMS 检索为效应矩阵

一次调用即可返回可直接分析的 (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：计算每个变异的预测器分数

从以下预测器来源中选择一种。选择仅会改变步骤 2；步骤 3–5 均相同。

预测器选项 A — ESM-C 6B SAE（完整示例）

对于每个变体，将残基窗口内的 SAE 激活值求和，计算 drop = max(0, WT − mut)，然后通过对 drop 取 top-K 均值，为每个变体汇总出一个分数：

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())

评估饱和突变扫描的两种方式——根据批次大小选择：

适用情况	使用方式	Forge 成本（19 个替换 × N 个位置）
饱和扫描包含 ≤100 个变体，或只需要每个变体的 top-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 成本减半；每次调用最多支持 100 个变体，足以覆盖单个位点的饱和诱变（19 个变体）或较短的位置扫描（例如位置 10-15：90 个变体）。对于更长的扫描，请拆分为多次批处理调用。

对于某些预测器分析所需的完整逐残基 × 逐特征张量，请改用循环模式。整个库的规模经过了充分测试：tests/integration/test_dms_pipeline_e2e_kras.py 会对 KRAS 的位置 10–25 运行约 300 个突变体。

预测器选项 B — AlphaMissense（hegelab 代理：分类结果，而非逐替换数值）

TU AlphaMissense 工具代理一个公共 API（hegelab.org），返回残基层面的分类结果，而非逐个（位置、替换氨基酸）的数值分数。提供三种调用方式；请根据规模选择：

工具	签名	返回值
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	一次调用即可获取完整蛋白质数据（用于全蛋白质 DMS 分析时成本最低）

这三种工具的响应结构相同。以 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 密钥。

预测器选项 C — 基于 ESM logits 的分数

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

预测器选项 C2 — ESM-2 masked-marginal（无需密钥，无需 API 密钥）

如果没有 ESM_API_KEY（选项 A/C 的工具需要该密钥），请使用 ESM2_score_missense_variant——它通过 HuggingFace 的免费 hf-inference 提供方运行 ESM-2，并针对一个错义变体返回 masked-marginal 对数似然比 logP(mut) − logP(wt)（Meier 等人，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 调用（没有批处理端点），因此对于饱和扫描，它比基于密钥的 ESM-C 批处理工具更慢——有密钥时优先使用选项 A/C；这是零配置的备用方案。长度超过约 1022 个残基的序列会自动围绕变体进行窗口截取（参见 metadata.windowed）。 LLR 是排序分数，而非经过校准的概率——步骤 3–5 会处理阈值判定，因此请将原始 LLR 直接输入 (20, n_positions) 矩阵。

预测器选项 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-retrieval skill 的元数据中获取 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 abundance 检测相反）。最简单的检查方法是计算 DMS 效应与已知方向固定的预测器之间的 Spearman 相关性（AlphaMissense 致病性分数始终为正值 = 破坏性更强）：

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} 执行检验，并报告全部三个结果——不同预测器的最佳 K 通常不同，因此应让比较过程透明，而不是针对结果调参。

MWU 是区分能力检验，但并非唯一有效的分析方法。 还可以考虑报告以下一项或两项：

• 预测器与 DMS 效应之间的 Spearman 秩相关性——可捕捉二元中性/破坏性划分所舍弃的分级校准信息。最适合回答“这个预测器的校准情况如何？”这类问题，而不只是“它能否区分？”。
• 在具有临床意义的破坏阈值处计算 AUROC（例如 ΔΔG > 1 kcal/mol）——可得出从 0 到 1 的数值，便于从业者理解，也可直接与文献中的预测器比较。

促使开发此 skill 的评估（KRAS folding 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）

---

比较两个预测器

在同一个 DMS 数据集上运行此技能两次（例如，SAE drops 和 AlphaMissense）， 使用相同的分层方式（disruptive_tail、neutral_abs、disruptive_quantile）， 然后进行比较：

指标	更好的预测器具备
p 值更低	区分能力更有把握
中位数差值更大（破坏性 − 中性）	效应量更大
覆盖范围中的 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 logits 评分	ESM_score_sequence
结构先验（用于预测器分析）	Structure_annotate_per_residue
下一步：针对每个热点分析机制	tooluniverse-residue-functional-mechanism-interpretation
最终可视化	tooluniverse-residue-functional-mechanism-interpretation 的第 7 步（带注释的热图 + 标注）