MetalKB：用团检测和统计势定位蛋白中的金属结合位点

赶了十多天软件，我回来了！

本文信息

标题：MetalKB：基于知识驱动图框架的蛋白金属结合位点预测
作者：Xuejun Zhao, Hao Li, and Sheng-You Huang*
发表时间：2026年3月25日（论文接收）
单位：华中科技大学物理学院，中国武汉
引用格式：Zhao, X., Li, H., & Huang, S.-Y. MetalKB: Predicting Metal Binding Sites on Proteins with a Knowledge-Based Graph Framework. Journal of Chemical Information and Modeling (2026). https://doi.org/10.1021/acs.jcim.6c00453
代码与资源：
- GitHub：https://github.com/huang-laboratory/MetalKB/
- 网页：http://huanglab.phys.hust.edu.cn/MetalKB/
- Zenodo：https://doi.org/10.5281/zenodo.18999183

摘要

金属离子在蛋白质的功能、调控和稳定性中发挥关键作用，因此，准确预测金属离子的结合位点，对于揭示相关生物过程的分子机制具有重要价值。本文提出了MetalKB，这是一种新的知识驱动框架，利用原子级统计势和图论策略来预测蛋白质上的金属离子结合位点。具体来说，先用clique检测算法识别可能的供体原子簇，并据此生成初始金属离子坐标；然后利用从蛋白—金属离子结合数据库推导得到的知识势，对这些候选坐标进行评估和局部细化；随后再通过空间距离阈值去除冗余预测。基于Metal3D和TEMSP提供的多样化基准数据集的评估表明，MetalKB在precision、recall和F1 score上与7种代表性方法相比具有有竞争力的表现，同时表现出较强的鲁棒性和参数稳定性。代表性结构案例进一步表明，MetalKB能够识别复杂的配位环境，包括多核金属位点和桥联金属位点。此外，它还能同时给出金属离子的三维坐标和残基级配位配体的预测。

abstract

图文摘要。这张图把MetalKB的主线压缩成了三个步骤：先从蛋白结构中提取候选供体原子，再用团检测找出可能共同配位的一组原子，最后用知识势打分与局部细化给出金属坐标和配位残基。

核心结论

MetalKB的核心创新，是把供体原子几何约束转写成图上的团检测问题，再用知识势进行筛选和局部优化。
本文既总结了Ca/Na、K、Mg和以Zn为代表的过渡金属几类主要配位偏好，又按各金属各自的数据集独立推导了金属特异性统计势。
在Metal3D与TEMSP两个基准上，MetalKB都给出了有竞争力的结果；第二篇会分别展开这两个基准各自关注的任务、评价标准和结果差异。
MetalKB的一个实际价值是同时给出金属离子的空间坐标和邻近配位残基，而不只是输出位点存在与否这一类粗粒度标签。

背景

金属离子在蛋白质中承担着多种角色，包括稳定结构、组织蛋白—蛋白界面、参与催化、调节信号转导以及维持离子稳态。已有研究估计，约30%–40%的蛋白需要一种或多种金属辅因子才能正常发挥功能，而锌尤其常见，在人体蛋白质组中约出现在10%的蛋白里。现有的金属结合位点预测方法大体可以分成序列方法和结构方法两类；结构方法内部又可以继续分成模板匹配、几何规则、机器学习和深度学习几条支线。

方法路线	代表方法	主要思路	当前进展	主要短板
基于序列的方法	MetalDetector、ZincFinder	从Cys/His富集片段、保守motif或序列特征出发预测金属结合位点	适合做大规模扫描，输入门槛低	缺少三维结构信息，较难描述空间上分散的配位残基
模板匹配结构方法	MIB、TEMSP	用已知金属位点模板、残基对模板或结构片段变换寻找相似配位环境	对已有模式通常表现较好，TEMSP在锌位点上尤其有代表性	依赖已知模板库，遇到新型配位模式或无同源模板时准确性容易下降
几何规则结构方法	BioMetAll、CHED	利用骨架预组织、几何triad、供体间距离和角度等规则筛位点	解释性强，容易与配位化学直觉对应	几何描述偏粗，难以完整表达金属—配体相互作用的能量差异
机器学习结构方法	ZincBindDB、PMM	从配位motif、结构特征和混合机器学习特征中学习位点模式，并常与几何过滤或打分策略结合	在特定金属或特定测试集上已经能取得较好表现，PMM是本文重点比较对象之一	特征工程和训练分布依赖明显，有些方法主要面向锌或过渡金属，跨金属泛化不稳定
深度学习结构方法	Metal3D、MetalSiteHunter、ESMBind、MasterOfMetals	分别利用体素化三维结构、三维卷积、蛋白预训练表示、图神经网络直接学习金属位点环境，并常与几何过滤或能量优化结合	代表了当前更自动化的数据驱动路线，Metal3D是本文重点比较对象；ESMBind、MasterOfMetals则体现了预训练模型和图神经网络在这一问题上的延伸	数据分布、金属类别覆盖和可解释性仍是主要瓶颈，对训练集组成较敏感
结构解析辅助方法	MIC	根据离子配位环境信息区分水分子和多种离子类型，面向冷冻电镜和晶体结构中的离子判别	更强调结构解释和离子归属，在实验结构解析场景中很有价值	关注点更偏离子类型判别与结构注释，不等同于从整条蛋白上完整搜索新的金属结合位点

实验上确定金属结合位点可以提供最直接的证据，但代价也高。质谱、X射线晶体学等技术可以提供高精度证据，不过成本高、周期长，不适合大规模筛选。因此，基于序列或结构的计算预测方法一直都很重要。问题在于，很多金属位点并不是线性序列上的连续motif，而是由空间上靠近、序列上相隔很远的残基共同构成，所以仅凭序列信息往往不足以描述真实配位环境。

结构方法虽然更接近真实配位环境，但也面临几个长期存在的问题：

模板法依赖已知模式，遇到新型配位环境或缺少合适模板的蛋白时，预测准确性就容易下降。
简单几何规则的信息量有限，距离和角度能描述一部分空间关系，却较难完整表达金属—配体相互作用。
QM/MM足够准但代价太高，不适合做常规的大规模扫描和筛选。

现有路线常见的问题是模板依赖过强、几何描述偏粗，或者训练分布过窄；一遇到多金属和复杂配位环境，泛化能力就容易受限。MetalKB用实验结构中已经积累的大量统计规律来做预测。

关键科学问题

怎样从整条蛋白结构里先找出合适的供体原子组合：真实金属位点通常至少包含3个配位供体，必须先把几何上可能同时配位的一组原子筛出来，否则后面的打分空间太大。
怎样把几何合理性与化学合理性结合起来：单靠供体—供体距离约束，可以筛掉很多明显不可能的情况，但仍会留下大量假阳性；还需要金属—原子相互作用势来进一步区分。
怎样兼顾多种金属类型而不过度依赖某一类训练集：Metal3D一类方法对锌表现突出，但对碱金属和碱土金属的泛化能力有限。MetalKB试图用金属特异性统计势缓解这个问题。
怎样处理多核和桥联位点：如果两个金属之间距离本来就很近，简单的空间聚类很容易把真实双核位点误删掉；方法必须能识别共享配体和近距离双金属构型。

创新点

把金属位点采样写成团检测问题，先用图论筛候选，再进入能量打分和细化。
从MESPEUS数据库推导距离依赖统计势，并与Lennard-Jones 12-6势混合，增强短程排斥和整体物理合理性。
显式引入羧酸侧链的虚拟供体节点，区分单齿、双齿、桥联等不同羧酸配位模式。
输出金属离子三维坐标与残基级配位信息，而不只是一个二分类标签。

MetalKB覆盖范围

这里要把测试覆盖范围和配位偏好的概括方式分开看。主文的多金属测试集明确包含Zn²⁺、Ca²⁺、Mg²⁺、Mn²⁺、Fe²⁺、Fe³⁺、Cu²⁺、Co²⁺、Ni²⁺、Na⁺和K⁺这11类金属离子；这些金属的统计势是按各自数据集独立推导的。只是从配位化学特征上看，它们又可以概括成4个代表类别：Ca/Na组、K组、Mg组，以及以Zn为代表的过渡金属组。Al³⁺、Mo、W这类离子没有出现在这篇的实际构建或测试范围里。

研究内容

fig1

图1：MetalKB的整体流程

阶段一：从金属蛋白结构中提取配位几何规则，并据此构建金属—蛋白原子对的知识势（knowledge-based potential / statistical potential）。
阶段二：先做基于clique的候选位点采样，再用混合势函数对候选位点评分、局部细化，并去除冗余预测。

MetalKB的整体思想是：先靠几何筛候选，再靠知识势与范德华势组成的混合势函数做化学判别。它比直接在整条蛋白上做均匀网格扫描更高效，因为大量非结合区域根本不会进入后续步骤。

基于clique的候选位点采样

fig2

图2：基于clique的候选位点采样

(a) 蛋白先被表示为供体原子集合，再转成图；节点是候选供体原子，只有当供体—供体距离落在统计得到的合理区间时，两点之间才连边。
(b)、(c)展示了羧酸氧参与金属配位时的两种典型模式，说明为什么仅靠均匀网格扫描较难区分这些模式。

这里的clique指的是完全连通子图。在MetalKB里，它表示一组供体原子两两之间都满足合理距离约束，因此有可能共同围成一个真实金属位点。

整个流程分成四步，而关键就在于先把搜索空间压缩到真正像配位簇的区域：

第一步，提取候选供体原子。过渡金属考虑Cys的SG、His的ND1/NE2、Glu的OE1/OE2、Asp的OD1/OD2；碱土和碱金属则考虑Asp、Glu、Asn、Gln、Ser、Thr的侧链氧以及所有残基的主链氧。
第二步，按供体—供体距离建图。对于过渡金属，两个供体原子距离落在2.4–5.2 Å时连边；其他金属类型则用图S1统计出来的各自区间，例如Ca²⁺和Mg²⁺是2.5–5.3 Å，K⁺是2.9–5.8 Å。
第三步，识别clique并做子团去冗余。这里要求clique至少包含3个供体原子；如果一个clique严格包含另一个较小clique，则保留大的那个，避免重复采样。
第四步，用供体原子几何质心作为初始金属坐标。这个坐标只作为后续局部精修的起点。

第四步里用到的几何质心写法是：

\[x = \frac{1}{n}\sum_{i=1}^{n} x_i,\qquad y = \frac{1}{n}\sum_{i=1}^{n} y_i,\qquad z = \frac{1}{n}\sum_{i=1}^{n} z_i\]

这里的$n$是clique中供体原子的数量，$(x_i, y_i, z_i)$是第$i$个供体原子的三维坐标。这个初始点是后续局部网格细化的起点：先识别出可能共同配位的一组供体，再用它们的几何中心给出初始金属位置。

候选位点的评分、局部细化与冗余去除

clique采样给出的几何质心只是一个初始候选位置，还不一定正好落在最低能量点上。为了进一步提高坐标精度，MetalKB又在每个初始坐标周围做了一轮局部网格细化：

初始坐标来自clique中供体原子的几何质心。
以这个初始坐标为中心，在2.5 Å半径内生成更密的局部候选点。
网格步长设为0.25 Å。
对每个候选点用前面定义的势函数逐一评分，保留能量最低的坐标作为最终预测位置。

这一轮细化的作用，是把初始几何估计进一步修到局部最低能量附近，从而改善金属坐标的定位精度。

去冗余时这里也特意避开了多核位点被误删的问题。多核金属簇里金属—金属距离多数在3–4 Å左右，因此MetalKB把冗余删除阈值设成2.5 Å。实际做法是先按能量从低到高排序，再检查预测点之间的距离；如果两个预测金属离子彼此小于2.5 Å，就保留能量更低的那个。另外，最终输出时只报告距离预测金属坐标4 Å以内的供体残基。例如锌位点只报告Cys、His、Glu、Asp这些符合统计规律的残基。

羧酸配位的特殊处理

Asp和Glu的羧酸基是这里最容易被低估的一类供体，因为一个残基上有两个氧原子，而这两个氧和金属的关系并不只有一种。羧酸可以形成四种不同的配位模式：

graph TB
    subgraph S1[单齿配位]
        O1((O<sub>1</sub>)) -->|配位键| M1((M<sup>n+</sup>))
        O2((O<sub>2</sub>))
    end

    subgraph S2[对称双齿配位]
        O3((O<sub>1</sub>)) -->|配位键| M2((M<sup>n+</sup>))
        O4((O<sub>2</sub>)) -->|配位键| M2
    end

    subgraph S3[非对称双齿配位]
        O5((O<sub>1</sub>)) -->|短键| M3((M<sup>n+</sup>))
        O6((O<sub>2</sub>)) -->|长键| M3
    end

    subgraph S4[桥联配位]
        O7((O<sub>1</sub>)) -->|配位键| M4a((M<sub>1</sub><sup>n+</sup>))
        O8((O<sub>2</sub>)) -->|配位键| M4b((M<sub>2</sub><sup>n+</sup>))
    end

    style M1 fill:#f9f,stroke:#333,stroke-width:2px
    style M2 fill:#f9f,stroke:#333,stroke-width:2px
    style M3 fill:#f9f,stroke:#333,stroke-width:2px
    style M4a fill:#f9f,stroke:#333,stroke-width:2px
    style M4b fill:#f9f,stroke:#333,stroke-width:2px

单齿配位：只有一个氧与金属形成配位键，另一个氧未参与。
对称双齿配位：两个氧同时与同一金属配位，键长基本相等。
非对称双齿配位：两个氧同时配位，但键长一短一长。
桥联配位：两个氧分别与两个不同的金属配位，使羧酸基成为连接两个金属中心的桥梁。

前三种模式都属于一个羧酸基围绕同一个金属中心配位的框架，只是两个氧参与的方式不同；桥联配位则更复杂，表示同一个羧酸基把配位关系延伸到两个金属中心。如果直接把两个羧酸氧都当成普通供体去做网格扫描，程序往往只能看到附近有两个氧，却分不清这里到底是一个氧单独配位，还是两个氧一起围住同一个金属，还是一个羧酸基同时参与两个金属中心。这也是普通均匀网格采样容易把这些几何模式混在一起的原因。

MetalKB的解决方案是在图论层面引入虚拟供体节点，专门代表把这个羧酸基当作一个双齿配位单元的情况：

graph TB
    subgraph G[羧酸基的图论表示]
        direction LR
        O1_node[氧原子O<sub>1</sub>]
        O2_node[氧原子O<sub>2</sub>]
        V_node[虚拟节点<br/>代表双齿配位]

        O1_node <-.无连接.-> O2_node
        O1_node <-.无连接.-> V_node
        O2_node <-.无连接.-> V_node

        O1_node -->|可与金属<br/>候选节点连边| Metal[金属候选位点]
        O2_node -->|可与金属<br/>候选节点连边| Metal
        V_node -->|可与金属<br/>候选节点连边| Metal
    end

    style V_node fill:#faa,stroke:#333,stroke-width:2px
    style O1_node fill:#afa,stroke:#333,stroke-width:2px
    style O2_node fill:#afa,stroke:#333,stroke-width:2px

这个虚拟节点是图模型里的占位符。关键在于两条不连边规则：

同一个羧酸基的两个氧之间不连边：它们不能同时出现在同一个clique中。
任意一个氧和它对应的虚拟节点之间也不连边：单齿模式与双齿模式互斥。

在clique搜索里，不连边就意味着这些节点不能同时出现在同一个clique里。因此，同一个羧酸基在一次候选采样中，只能以一种方式被表示：

选其中一个氧 → 当作单齿配位的一侧供体。
选虚拟节点 → 当作双齿配位单元。
不会把两个氧和虚拟双齿节点同时塞进同一个候选簇里，避免重复计数。

这套表示方式首先解决的是单齿与双齿配位的区分问题。桥联配位的关键不是两个氧一定同时进入同一个clique，而是同一个羧酸基可以被不同金属中心共享。虚拟节点和互斥规则有助于避免把同一个羧酸基在单个候选簇中重复计数；多核位点能否被保留，则主要取决于后续2.5 Å冗余阈值没有把彼此真实相邻的金属中心误合并。

这样设计的核心是在图论采样阶段，先把羧酸这种多模式供体拆成几种互斥的几何表示。这样后面的clique搜索和势函数打分才能分清：这是单齿、双齿，还是可能出现在更复杂的桥联环境里。

这里还有一个正文没有写透、但源码补上的实现细节。论文明确说，虚拟供体节点是加在两个羧酸氧之间，用来表示潜在双齿配位；SI也补充说明，程序里因为引入了这种virtual-atom representation，所以同一羧酸基两个氧之间的距离不纳入donor-donor距离统计。

源码进一步表明，这个虚拟节点的坐标就是两个羧酸氧坐标的几何中点。
同时，它没有被赋予单独的metal-virtual atom势函数类型，主要用于clique采样阶段的几何表示与质心估计。
这个设计把“羧酸基按双齿单元参与配位”的几何信息带进候选簇；后续进入金属—原子打分和残基输出的，仍然是实际存在的蛋白原子。

知识势与混合势函数

本文里的$w_{ij}(r)$是从结构数据库统计出来的知识势，英文常写作knowledge-based potential，也就是一类statistical potential；真正用于打分的是再与范德华势组合后的混合势函数$u_{ij}(r)$。知识势$w_{ij}(r)$基于观测到的原子对距离分布：某种金属—原子相互作用在实验结构中出现得越频繁，对应的能量就越低。

MESPEUS是这里的主要数据来源。这个数据库专门整理蛋白中的金属位点，且只收录分辨率优于2.5 Å、由X射线晶体学或冷冻电镜解析的结构，不包含NMR或分辨率不明的条目，因此适合做几何统计。文中先统计不同金属偏好的供体类型。Table 1给出的是残基供体频率图谱。为了让这个统计更直观，可以把正文里的信息压缩成下面这张总结表：

金属类别	主要高频供体	文章强调的配位特征
Ca²⁺、Mg²⁺	Asp/Glu羧酸氧、主链羰基氧	偏好氧供体，Mg²⁺的配位特征与Ca²⁺仍有明显差别
Na⁺、K⁺	主链羰基氧和各类侧链氧	以氧供体为主，K⁺对主链羰基氧尤其常见
Mn、Fe、Co	His咪唑氮、Asp/Glu羧酸氧	兼具His偏好与对酸性残基的明显使用
Ni、Cu、Zn	His咪唑氮、Cys硫原子、Asp/Glu羧酸氧	His/Cys偏好最突出，尤其Cu、Zn对Cys的偏好很明显

真正影响后续建模的是，文中先根据供体组成、离子半径和配位特征的共性，把主要配位偏好概括成Ca/Na组、K组、Mg组，以及以Zn为代表的过渡金属组；但在具体实现里，统计势仍然按各金属各自的数据集独立推导，并没有在不同金属之间共享参数。
为了降低冗余，这里还用CD-HIT在30%序列一致性阈值上做了去冗余。最终用于知识势推导的数据量分别是：Zn结合蛋白2568个、Ca2375个、Mg3451个、K778个，这意味着这些势函数建立在规模较大的结构统计样本之上。

知识势函数（Eq 2）

本文使用逆Boltzmann形式把局部富集程度转换成相对能量：

\[w_{ij}(r) = -k_B T \log \left[ \frac{\rho_{ij}^{\mathrm{obs}}(r)}{\rho_{ij,\mathrm{bulk}}^{\mathrm{obs}}} \right]\]

这里，$\rho_{ij}^{\mathrm{obs}}(r)$是金属离子$i$与原子类型$j$在距离$r$处的观测数密度，$\rho_{ij,\mathrm{bulk}}^{\mathrm{obs}}$是参考球体中的平均背景数密度，可以把它理解成没有特殊配位偏好时的基线水平。

这里的感觉和RDF很接近：本质上都是在比较某个距离壳层里的局部富集程度，只不过本文把参考背景写成了10 Å参考球内的平均数密度。计算时把$k_B T$设为1，所以这个式子更接近一种相对打分势，重点是比较不同相互作用在结构数据库里出现得是否异常频繁。

混合势函数（Eq 5）

本文把知识势和范德华势拼在一起：

\[u_{ij}(r)= \begin{cases} \min \left[w_{ij}(r),\, v_{ij}(r)\right], & r \le 3.0\ Å \\ \dfrac{v_{ij}(r)e^{-v_{ij}(r)} + w_{ij}(r)e^{-w_{ij}(r)}}{e^{-v_{ij}(r)} + e^{-w_{ij}(r)}}, & r > 3.0\ Å \end{cases}\]

这里的$v_{ij}(r)$是Lennard-Jones 12-6势。这个分段形式的意义很明确：在3.0 Å以内，直接取知识势和范德华势里更保守的那个，避免短程碰撞被知识势低估；在3.0 Å以外，再用指数加权把两者平滑拼接起来，让长程能量自然衰减到0。

由于权重写成了$e^{-v_{ij}(r)}$和$e^{-w_{ij}(r)}$，所以能量更低的那一项会占更大权重。这更接近一种低能优先的平滑拼接，既保留了实验结构统计里的配位偏好，又不会在短距离给出明显不合理的能量形状。

总打分时，会对金属离子与周围相关蛋白原子对的相互作用逐一求和。真正参与打分的是混合势函数$u_{ij}(r)$：它由知识势$w_{ij}(r)$与Lennard-Jones 12-6势$v_{ij}(r)$组合得到，在短距离保留保守的排斥与势阱形状，在较长距离则通过低能项权重更高的指数加权实现平滑过渡。这个设计的重点是同时保留配位偏好与短程排斥。

fig3

图3：四种代表性原子对势函数。图3把前面的统计规律落实到了具体势函数上，也让混合势函数的形状是否合理变得更直观。

(a) Zn−S.3的势阱最深、最窄，说明Zn与Cys硫原子的配位更强、更刚性，这与它常承担结构稳定作用一致。
(b) Zn−N.ar在约2.0 Å和4.1 Å处出现两个势阱，前者对应组氨酸咪唑氮与Zn的直接配位，后者则更浅，反映出更远距离下的立体效应。
(c) Zn−O.co2的势阱较宽，可从约2.0 Å延伸到4.0 Å，体现了羧酸氧既可以单齿配位，也可能通过双齿或桥联方式参与配位。
(d) Ca−O.co2在约2.3 Å附近达到极小值，并在较远距离保留次级势阱，说明Ca与羧酸氧的配位几何和Zn并不相同。

深入思考：

知识势和Lennard-Jones势的来源并不相同，这里的目标也不是得到严格可比的绝对相互作用能，而是构造一个用于排序和定位的有效势函数。知识势来自逆Boltzmann关系，反映实验结构里某种金属—原子相互作用在不同距离上的相对偏好；范德华项补上了短程排斥和势阱形状的物理约束。到了更远距离，知识势会随着局部密度接近背景而趋近0，Lennard-Jones项也会快速衰减到0，因此长程行为是一致的。

原文并没有单独给出一张“知识势和范德华势数量级逐点对齐”的标定图，Figure 3提供的是一种间接验证：如果拼接后的势函数在极小值位置、势阱深浅和曲线宽窄上都符合常见配位化学规律，那么这个混合方式至少在方法上是自洽的。Figure 3展示了“这个混合势函数像不像一个合理的配位势”。从这四条曲线看，最小值位置与常见配位键长基本一致，短程没有出现明显不合理的深井，长程也能逐步回到0。这就是本文采用这种拼接方式的主要依据。

Mendelevium

Contact