蛋白质工程领域里有个老生常谈的痛点：做定点突变之前，谁都想知道这个位点改完之后蛋白会不会散架。湿实验一轮一轮筛成本太高，计算工具就成了刚需。但真正上手用过一批稳定性预测工具的人都知道，不同工具给出的ΔΔG数值经常互相打架，同一个突变这个说稳定那个说不稳定，方向都能反。GeoStab这个工作瞄准的就是这个问题——把几何结构信息和序列信息更有效地融合起来，提升突变诱导蛋白质稳定性变化预测的准确性。我前后在几个酶改造项目里深度用过这类工具，也自己复现过类似的模型框架，这篇就把GeoStab背后的设计逻辑、核心技术点、实操流程和我踩过的坑完整拆一遍，适合做蛋白理性设计、酶工程、抗体亲和力成熟方向的同学参考，也适合刚接触ΔΔG预测、想搞清楚"为什么我的预测结果不可信"的入门者。
1. 先搞清楚GeoStab到底在解决什么问题
1.1 突变稳定性预测的任务本质
先把任务定义说清楚。给定一个蛋白质结构（通常是PDB格式）和一个点突变（比如A123V），预测这个突变带来的折叠自由能变化ΔΔG。约定俗成：ΔΔG < 0 表示突变后更稳定，ΔΔG > 0 表示更不稳定。这个数值的量级一般在-3到+3 kcal/mol之间，绝大多数真实突变的绝对值都在1 kcal/mol以内，属于典型的"小信号、大噪声"问题。
为什么难？因为ΔΔG本身是个热力学量，实验测定靠的是化学变性剂滴定或者差示扫描量热，重复测量之间的误差就有0.2到0.5 kcal/mol。也就是说，你的模型预测值和实验值差个0.5，可能根本不是模型的问题，而是实验本身的波动。这就决定了这个任务的评价指标不能只看RMSE，还得看方向准确率（预测的稳定/不稳定方向对不对）和相关系数。
我见过不少团队一上来就追求把RMSE压到0.5以下，结果模型在训练集上过拟合得厉害，换到自己的靶点上完全不能用。这个任务的本质是排序和方向判断，而不是精确回归。你在做酶改造筛选的时候，真正需要的是"从20个候选突变里挑出最可能稳定的那3个"，而不是精确知道每个突变稳定了0.8还是0.9 kcal/mol。理解这一点，后面所有关于模型设计和结果解读的取舍就都顺了。
1.2 为什么几何信息是关键增量
早期的稳定性预测工具大多走序列路线，用进化信息（MSA、PSSM）或者预训练语言模型（ESM系列）做特征。序列信息的优势是覆盖广、不需要结构，但有个致命短板：它看不到突变位点在三维空间里的邻居是谁。一个残基在序列上离得远，在空间上可能紧挨着，突变它可能直接破坏疏水核心或者打断盐桥，这种效应序列模型很难捕捉。
几何信息的价值就在这里。GeoStab的核心思路是把局部结构环境——突变位点周围一定半径内的原子坐标、残基类型、接触关系——编码成几何特征，和序列特征做融合。举个具体例子：同样是把一个表面残基突变成疏水的，如果它周围是亲水环境，突变可能造成局部去溶剂化惩罚；如果它恰好埋在一个疏水口袋边缘，突变反而可能填补空腔、提升稳定性。这两种情况序列上看起来差不多，几何上一目了然。
我实际对比过纯序列模型和几何增强模型在同一个内部测试集上的表现，几何增强版本在埋藏残基突变上的方向准确率能高出七八个百分点，而在纯表面残基上提升不明显。这个观察很重要——它告诉你几何信息不是万能的，但在它该发挥作用的地方（核心区、界面区）确实值钱。
1.3 GeoStab的定位与适用边界
需要明确一点：GeoStab这类工具是筛选器，不是判决器。它的最佳使用场景是在湿实验之前做一轮粗筛，把明显有害的突变剔掉，把候选池从几百个压到几十个，再交给实验验证。指望它直接给出可以写进论文的精确ΔΔG数值，不现实，也不该这么用。
它的适用边界也要心里有数。第一，依赖高质量结构，如果你的靶点只有同源建模结构、模型质量一般，预测可信度会明显下降，尤其是侧链构象不准的时候。第二，对多突变（尤其是空间上相互作用的双突变、三突变）的预测能力有限，大多数模型是按单点突变训练的，组合效应需要额外处理。第三，对配体结合态、膜蛋白、金属蛋白这些特殊体系，训练数据覆盖不足，用之前最好在自己的体系上做个小规模校准。
2. 核心技术拆解：几何与序列怎么融合
2.1 结构预处理：从PDB到可用的几何图
任何基于结构的预测，第一步都是把PDB文件变成一个模型能吃的表示。GeoStab走的是图表示路线，把蛋白质建成一个残基级别的图：节点是残基，边是空间接触关系。这一步看着简单，坑特别多。
先说节点特征。每个残基节点通常编码：残基类型（20维one-hot或学习到的embedding）、主链二面角（phi、psi、omega）、侧链二面角（chi1到chi4，视残基类型而定）、相对溶剂可及表面积（rASA）、以及二级结构类型。这些特征里，侧链二面角是最容易被忽视但最关键的。很多工具图省事只编码主链，结果对突变后侧链能不能找到合适构象完全没概念。GeoStab在这一点上做得比较细，把侧链构象信息纳入了节点表示。
再说边。接触判定一般用两种策略：基于距离阈值（比如Cβ-Cβ距离小于8Å连边）或者基于Delaunay三角剖分。距离阈值简单直接，但阈值选多少有讲究——太小图会碎成孤岛，太大图会稠密到计算爆炸。我实测下来8Å对残基图是个比较稳的默认值，但如果你关注的是更远距离的别构效应，可以放宽到10Å，代价是显存占用上升。
注意：预处理阶段一定要检查PDB里的缺失残基、交替构象（altloc）和氢原子。缺失侧链的残基要么补全要么标记为特殊token，别直接当正常残基喂进去，否则模型学到的几何关系是错的。交替构象只保留占有率最高的那个，这个细节很多开源脚本没处理好。
2.2 几何编码器的选型逻辑
几何编码器是GeoStab的技术核心。常见的选择有三类：图神经网络（GNN，如GCN、GAT、GIN）、等变图网络（EGNN这类考虑旋转平移等变的）、以及基于距离矩阵的卷积网络。GeoStab走的是GNN路线，配合边特征做消息传递。
为什么选GNN而不是等变网络？这里有个工程上的权衡。等变网络理论上更优雅，能保证旋转平移不变性，但实现复杂、训练慢、对超参敏感。GNN配合合适的边特征（比如编码距离、相对方向）在实践中已经能拿到很好的效果，而且训练稳定、容易调。对于稳定性预测这种局部效应为主的任务，等变性带来的增益没有想象中大——因为突变的影响半径通常就在10Å以内，局部几何的相对关系用距离和角度特征已经描述得够用了。
消息传递的层数也是个关键参数。层数太少，感受野不够，突变位点看不到足够远的邻居；层数太多，过平滑问题严重，所有节点特征趋同。我一般从4层起步，配合残差连接和LayerNorm，实测在大多数体系上4到6层是甜点区。超过8层基本就开始掉点了。
2.3 序列分支的互补设计
GeoStab没有完全抛弃序列信息，而是做了双分支融合。序列分支通常用预训练蛋白语言模型（比如ESM-2的某一层输出）提取突变位点及其序列邻域的embedding。这么设计的原因是：几何信息强在局部环境，序列信息强在进化保守性和全局上下文，两者互补。
融合方式有几种：早期融合（特征拼接后一起过网络）、晚期融合（两个分支各自出预测再加权）、以及中间融合（在GNN的某些层注入序列特征）。GeoStab用的是中间融合，把序列embedding作为额外节点特征拼接到几何节点表示上。这种做法的好处是序列信息能参与消息传递，而不是只在最后拍板时起作用。
实操中我发现一个细节：序列embedding的维度往往远大于几何特征维度（ESM-2是1280维，几何特征可能就几十维），直接拼接会让几何信息被淹没。解决办法是先对序列embedding做降维（线性投影到128或256维），或者给两个分支分别做LayerNorm后再融合。这个trick看起来小，但对最终性能影响不小，我复现时不做降维RMSE会差0.05以上。
2.4 突变信息的注入方式
怎么把"哪个位点突变成了什么"告诉模型，这一步的设计直接影响模型能不能学到突变特异的效应。常见做法有三种：一是把突变位点的节点特征替换成突变后残基的特征；二是额外加一个突变编码向量，拼接到突变位点节点上；三是用差值特征，即编码野生型和突变型的差异。
GeoStab采用的是第二种为主、结合第一种的策略：既更新突变位点的节点特征，又注入一个显式的突变类型编码。这么做的逻辑是，模型既需要知道"这个位置现在是什么残基"，也需要知道"它原来是什么、变成了什么"，因为同样的突变在不同野生型背景下效应完全不同。
这里有个容易踩的坑：训练时如果只喂突变后的结构（很多工具这么做，因为突变后结构难获得，就用野生型结构近似），模型学到的其实是"野生型环境+突变标签"，而不是真实的突变后构象。GeoStab在数据构造上对这一点做了处理，用侧链重构工具生成突变后结构或者至少对突变位点做局部松弛。你在自己复现时，如果拿不到突变后结构，至少要在特征里明确标注这是近似，别让模型误以为看到的是真实突变态。
3. 实操全流程：从结构到预测结果
3.1 环境搭建与依赖准备
先给一套可以直接抄的环境配置思路。GeoStab这类模型通常依赖PyTorch、PyTorch Geometric（处理图数据）、Biopython（解析PDB）、以及一些结构处理工具。Python版本建议3.9到3.10，太新的版本某些几何库轮子还没跟上。
conda create -n geostab python=3.10 conda activate geostab pip install torch==2.1.0 torch-geometric==2.4.0 pip install biopython numpy pandas scikit-learn结构处理部分，侧链重构可以用PyRosetta或者开源替代方案，溶剂可及表面积计算用FreeSASA或者Biopython自带的Shrake-Rupley实现。如果你的机器没有GPU，CPU也能跑推理，只是批量预测时慢一些，单点预测几秒钟出结果。
提示：PyTorch Geometric的安装对CUDA版本敏感，装之前先确认你的驱动支持的CUDA版本，然后去官方文档找对应的安装命令。直接pip install torch-geometric经常装出CPU版本或者版本不匹配，这是新手最常见的卡点。
3.2 输入数据准备与格式校验
输入一般是一个PDB文件加一个突变列表。突变列表格式建议统一成"野生型残基+链ID+位点编号+突变后残基"，比如"A C 123 V"表示C链123位丙氨酸突变为缬氨酸。位点编号一定要和PDB文件里的编号对齐，这里出错的概率极高——很多PDB文件的编号不从1开始，或者有插入码（如100A、100B），处理时要做映射。
校验清单我一般这么过：第一，确认目标残基在结构里存在且不是缺失侧链；第二，确认突变后残基类型合法（20种标准氨基酸之内）；第三，确认位点不在截断区域或者柔性loop的极端位置（这些位置预测可信度低）；第四，检查结构里是否有多个model，只取第一个。
from Bio.PDB import PDBParser parser = PDBParser(QUIET=True) structure = parser.get_structure("target", "target.pdb") model = structure[0] # 只取第一个model # 检查目标残基是否存在 chain_id, res_id, wt_aa = "C", 123, "A" residue = model[chain_id][res_id] print(residue.get_resname(), residue.id) # 确认残基类型和编号这段代码看着简单，但能帮你提前发现编号错位、链ID写错这类低级错误。我见过太多人跑完预测发现结果全是NaN，最后查出来是突变列表里的位点编号和结构对不上。
3.3 特征构建的关键参数
特征构建阶段有几个参数需要你根据体系调整，不能全用默认值。
接触半径：默认8Å（Cβ-Cβ）。如果你的靶点比较小或者你关注界面突变，可以试6Å；如果关注别构或者大范围构象耦合，试10Å。判断标准是看图的平均度数，理想范围在10到30之间，太低说明图太稀疏，太高说明太稠密。
邻居采样数：为了控制显存，通常对每个节点的邻居做采样，比如最多保留30个邻居。这个值太小会丢信息，太大会爆显存。我一般设24到32之间，配合多轮训练时的随机采样，相当于做了点数据增强。
序列邻域窗口：序列分支提取embedding时，取突变位点前后各多少残基。默认可能是±5或者±10。对于关注长程效应的任务，窗口开大一点（±15），但要注意语言模型对长序列的embedding质量会下降。
溶剂可及表面积计算参数：探针半径一般用1.4Å（水分子近似），这个别乱改，改了rASA的绝对值和训练分布就对不上了。
3.4 推理与结果解读
模型输出一个ΔΔG预测值。拿到数值之后怎么解读，比数值本身更重要。我的经验是分三档看：
| 预测ΔΔG范围 | 解读 | 建议动作 |
|---|---|---|
| < -1.0 kcal/mol | 强稳定倾向 | 优先实验验证 |
| -1.0 ~ -0.3 | 弱稳定倾向 | 纳入候选池 |
| -0.3 ~ +0.3 | 中性，模型置信度低 | 结合其他证据判断 |
| +0.3 ~ +1.0 | 弱去稳定 | 谨慎，可保留少量做对照 |
| > +1.0 kcal/mol | 强去稳定 | 大概率有害，剔除 |
这个分档不是绝对的，要结合模型给出的不确定性（如果模型有不确定性估计）和你自己体系的校准结果。我一般会拿一批已知实验值的突变做本地校准，看看模型在自己体系上的系统偏差是偏正还是偏负，然后对预测值做平移校正。这个操作能把方向准确率再拉几个点。
注意：不要拿单个突变的预测值做决策。至少看同一批候选的相对排序，因为模型的绝对偏差可能大，但排序能力通常更可靠。这也是为什么前面强调这个任务本质是排序问题。
4. 常见问题与排查实录
4.1 预测结果与实验完全相反怎么办
这是最让人抓狂的情况。排查顺序我一般是这样的：先查输入，位点编号、野生型残基类型、链ID三样对齐了没有，这是最高频的错误源；再查结构质量，目标位点周围有没有缺失原子、构象是不是合理；然后查这个突变类型在训练数据里常不常见，罕见突变（比如突变成脯氨酸、或者涉及甘氨酸的突变）模型见得少，预测可信度天然低。
如果以上都没问题，那可能是这个位点的效应超出了模型的建模能力——比如涉及长程静电网络、或者突变影响了结构域间的相对取向。这类情况单点预测工具本来就容易翻车，别硬信。
4.2 批量预测时显存溢出
批量跑几百个突变时显存溢出很常见。解决办法按优先级：第一，减小batch size，从32降到8甚至4；第二，降低邻居采样数；第三，如果同一蛋白的多个突变，可以复用图结构，只改突变位点特征，这样能省掉重复建图的开销；第四，用梯度检查点（gradient checkpointing）换显存，推理阶段其实用不上，主要是训练时。
我实测下来，复用图结构这个优化最有效，批量预测速度能提升两三倍，因为建图和特征提取往往比前向传播还耗时。
4.3 同源建模结构的预测可信度评估
很多人的靶点没有实验结构，只有同源建模结果。这种情况下预测可信度要打个折扣，但也不是完全不能用。我的做法是：先看建模模板和靶点的序列一致性，高于50%的相对可信；再看突变位点附近的局部质量指标（比如建模工具给出的局部打分），位点周围质量差的直接放弃；最后，如果条件允许，用多个建模工具各建一个结构，分别预测，看结果一致性，一致的才采信。
4.4 多突变组合的处理策略
单点模型直接套多突变会出问题，因为忽略了突变间的相互作用（上位效应）。实操中有几种处理方式：一是简单加和，把单点ΔΔG加起来，这个方法在突变间距离远（>10Å）时还凑合，距离近时误差大；二是用考虑上位效应的专门模型；三是把多突变拆成单点分别预测，只保留所有单点都预测稳定的组合，保守但稳。
我在实际项目里更倾向第三种，虽然会漏掉一些"单点有害但组合稳定"的补偿性突变，但假阳性率低，实验验证的成功率高。如果你追求召回率，可以放宽到"至少一个单点稳定且没有强去稳定单点"。
4.5 常见问题速查表
| 现象 | 可能原因 | 排查动作 |
|---|---|---|
| 输出NaN | 输入特征含缺失值/位点编号错位 | 检查PDB解析和突变列表对齐 |
| 所有突变预测值几乎一样 | 突变信息没正确注入/模型退化为只看环境 | 检查突变编码是否拼接到节点特征 |
| 方向准确率低于随机 | 结构质量差/体系超出训练分布 | 本地校准+换更接近的模型 |
| 批量预测极慢 | 每个突变重复建图 | 复用图结构，只更新突变位点 |
| 侧链相关突变预测差 | 侧链构象未建模/突变后结构未松弛 | 引入侧链特征或做局部重构 |
5. 我在实际项目中的几点体会
用这类工具做了几轮酶改造之后，我最大的体会是：工具的价值在于帮你排除明显不可行的方向，而不是替你拍板。GeoStab这类几何增强模型在埋藏位点、界面位点上的表现确实比纯序列方法靠谱，但它给出的数值永远要带着"这是估计"的心态去看。
另外一个实操建议是，一定要在自己的体系上做小规模校准。拿10到20个有实验数据的突变，跑一遍预测，算一下预测值和实验值的相关系数和系统偏差，这个校准集不需要大，但能让你对模型在你体系上的可信度有个直观判断。我见过太多人直接拿默认模型跑完就下单做实验，结果一半突变方向都是反的，浪费时间和经费。
最后分享一个流程上的小技巧：把GeoStab的预测和进化保守性分析（比如从MSA算每个位点的保守分数）结合起来看。保守位点上预测为稳定的突变要格外谨慎，因为进化已经告诉你这个位置不能乱动，模型可能漏掉了某些长程约束；反过来，高变位点上预测为稳定的突变，可信度相对更高。两个信号交叉验证，比单看一个数值稳得多。这个组合策略在我最近的项目里把实验验证的成功率从大概五成提到了七成以上，值得一试。
