蛋白质工程领域里有个老生常谈的痛点做定点突变之前谁都想知道这个位点改完之后蛋白会不会散架。湿实验一轮一轮筛成本太高计算工具就成了刚需。但真正上手用过一批稳定性预测工具的人都知道不同工具给出的ΔΔ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 python3.10 conda activate geostab pip install torch2.1.0 torch-geometric2.4.0 pip install biopython numpy pandas scikit-learn结构处理部分侧链重构可以用PyRosetta或者开源替代方案溶剂可及表面积计算用FreeSASA或者Biopython自带的Shrake-Rupley实现。如果你的机器没有GPUCPU也能跑推理只是批量预测时慢一些单点预测几秒钟出结果。提示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(QUIETTrue) 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算每个位点的保守分数结合起来看。保守位点上预测为稳定的突变要格外谨慎因为进化已经告诉你这个位置不能乱动模型可能漏掉了某些长程约束反过来高变位点上预测为稳定的突变可信度相对更高。两个信号交叉验证比单看一个数值稳得多。这个组合策略在我最近的项目里把实验验证的成功率从大概五成提到了七成以上值得一试。