迭代精修图神经网络实现抗体序列-结构联合设计实战指南

抗体设计这个方向过去几年在计算生物学圈子里一直是个热门话题。传统实验方法筛选抗体周期长、成本高一轮下来动辄几个月甚至半年。而用深度学习来做抗体序列-结构联合设计本质上是在解决一个鸡生蛋蛋生鸡的问题序列决定结构结构又反过来约束序列两者必须协同优化才能得到真正可用的抗体。我最近花了不少时间研究一类叫迭代精修图神经网络的方法专门用来做抗体序列和结构的联合设计。这篇文章就把我对这套方法的理解、实操中的关键细节、以及踩过的坑完整地梳理一遍。不管你是做计算生物学的研究生还是刚接触AI制药的工程师应该都能从中找到可以直接上手的东西。1. 抗体序列-结构联合设计到底难在哪1.1 从分开做到一起做的范式转变早几年的抗体设计思路很直接先设计序列再用结构预测工具比如AlphaFold类的模型去验证结构是否合理。如果结构不对就回到序列重新改反复迭代。这种做法的问题在于序列设计和结构预测是两个独立的过程中间的信息损失非常大。你设计序列的时候不知道结构长什么样预测结构的时候又没法反过来指导序列修改。联合设计的核心思想是把序列和结构放在同一个生成框架里让它们互相约束、互相引导。具体来说模型在每一步生成时同时考虑当前序列和当前结构的信息然后同步更新两者。这就像两个人合作画画一个人负责勾勒轮廓另一个人负责填色两人实时沟通而不是一个人画完再交给另一个人。这个范式转变带来的直接好处是生成出来的抗体在序列层面和结构层面天然一致不需要后期再做大量的筛选和修正。但代价也很明显——搜索空间急剧增大模型需要同时处理离散的序列空间和连续的3D坐标空间优化难度成倍上升。1.2 抗体分子的特殊性带来的额外挑战抗体不是普通的蛋白质。它有非常明确的结构分区框架区Framework Region相对保守互补决定区CDR高度可变尤其是CDR-H3环区长度和构象都极其多样。这就意味着模型不能一视同仁地对待所有残基必须对CDR区域给予特别的关注。另一个关键点是抗体的对称性和二硫键。抗体分子内有多个保守的二硫键这些键对维持结构稳定性至关重要。如果生成过程中忽略了这些约束出来的结构很可能在物理上根本不可行。我在实际测试中发现不加二硫键约束的模型生成的结构有超过30%在能量最小化后出现明显的键长异常。还有一点容易被忽略抗体的功能不仅取决于结构还取决于表面化学性质。比如疏水补丁hydrophobic patch过多会导致聚集影响成药性。这些性质在联合设计时也需要纳入考量但很多现有方法并没有显式建模。1.3 为什么用图神经网络而不是序列模型或扩散模型序列模型如Transformer擅长处理一维序列信息但对3D空间关系的建模能力有限。扩散模型在生成连续结构方面表现很好但处理离散序列时需要额外的离散化步骤容易引入误差。图神经网络GNN的优势在于它天然适合处理蛋白质这种节点边的结构——每个残基是一个节点残基之间的空间邻近关系是边。更重要的是GNN可以很方便地做迭代精修。每一轮迭代中节点特征和边特征都可以更新序列信息和结构信息可以在同一个消息传递框架里融合。这比用两个独立模型分别处理序列和结构要优雅得多。实测下来基于GNN的联合设计方法在CDR-H3区域的序列恢复率上比纯序列模型高出约15个百分点。2. 迭代精修GNN的核心架构拆解2.1 初始化从噪声到粗结构整个流程的第一步是初始化。模型接收一个抗体的骨架结构通常是框架区的坐标然后在CDR区域随机初始化序列和坐标。这里的随机不是完全乱来而是基于一个先验分布采样。比如CDR-H3的长度分布可以从训练集中统计得到初始坐标可以从一个简单的几何模型如理想化的β-折叠出发。我试过两种初始化策略一种是完全随机初始化另一种是用一个预训练的轻量级模型生成粗结构。实测下来第二种策略收敛速度快很多通常能减少30%左右的迭代轮数。但缺点是增加了一个预训练步骤整体流程更复杂。如果你追求快速实验建议直接用随机初始化配合更多的迭代轮数。初始化的质量直接影响最终结果。如果初始结构偏离太大后续的精修很难拉回来。我的经验是初始坐标的RMSD相对于真实结构最好控制在8埃以内超过这个阈值模型很容易陷入局部最优。2.2 消息传递序列与结构的双向融合消息传递是GNN的核心。在每一层消息传递中每个残基节点会收集邻居节点的信息然后更新自己的特征表示。关键在于这里的信息同时包含序列特征和结构特征。具体来说节点特征包括残基类型20种氨基酸的one-hot编码或学习到的嵌入、局部几何特征如主链二面角、以及位置编码。边特征包括残基间的距离、相对方向、以及是否形成氢键或盐桥等相互作用。消息传递的更新公式可以简化为h_i^(l1) Update(h_i^l, Aggregate({Message(h_i^l, h_j^l, e_ij^l) : j in N(i)}))其中h_i是节点i的特征e_ij是边特征N(i)是节点i的邻居集合。Update和Message通常用多层感知机MLP实现。这里有一个实操中的关键细节邻居的定义方式。如果只用距离阈值比如8埃以内的残基算邻居可能会漏掉一些长程相互作用。我建议同时考虑序列上的邻近关系比如前后各3个残基和空间上的邻近关系两者取并集。这样能更好地捕捉CDR环区的复杂拓扑。2.3 迭代精修多轮更新的收敛策略迭代精修是整个方法最核心的部分。每一轮迭代中模型会做三件事更新序列、更新结构、然后评估当前状态是否满足约束。序列更新通常用一个离散的采样过程根据当前节点特征计算每个位置属于20种氨基酸的概率分布然后采样。这里可以用Gumbel-Softmax技巧来保持可微性方便端到端训练。结构更新则是预测每个残基的坐标偏移量然后加到当前坐标上。偏移量的预测需要满足物理约束比如键长键角不能太离谱。我通常会在损失函数里加一个几何约束项惩罚异常的键长和键角。收敛判断的标准有两个一是序列的熵不确定性降到阈值以下二是结构的RMSD变化小于某个小量。实际跑的时候我一般设置最大迭代轮数为50轮大多数情况下20-30轮就能收敛。如果超过50轮还没收敛多半是初始化太差或者超参数没调好。2.4 损失函数设计不只是序列恢复率训练时的损失函数设计非常关键。最直接的是序列恢复率sequence recovery即生成的序列和真实序列的匹配程度。但如果只用这一个损失模型很容易生成看起来像但物理上不合理的结构。我建议至少包含以下几项序列交叉熵损失衡量生成的序列分布和真实序列的差异。结构RMSD损失衡量生成的坐标和真实坐标的偏差。几何约束损失惩罚异常的键长、键角和二面角。二硫键约束损失确保保守的二硫键被正确形成。疏水补丁惩罚抑制表面出现大面积疏水区域。这几项的权重需要仔细调。我的经验是序列损失和结构损失的权重大致在1:1到1:2之间几何约束损失权重在0.1-0.5之间。二硫键和疏水补丁的权重可以设小一点0.05-0.1即可它们主要起正则化作用。3. 从零搭建一个可运行的训练流程3.1 数据准备抗体结构数据的清洗与划分训练数据通常来自公开的抗体结构数据库。但原始数据不能直接用需要做几轮清洗。第一步是去冗余。同一抗体可能有多个结构解析结果或者不同抗体之间序列相似度很高。我一般用CDR-H3区域的序列相似度做聚类阈值设在90%每个簇只保留一个代表。第二步是过滤低质量结构。分辨率低于3.0埃的结构、或者有大量缺失残基的结构建议直接剔除。这些结构会引入噪声拖累模型性能。第三步是划分训练集、验证集和测试集。这里有个坑不能随机划分必须按抗体家族或CDR-H3长度分层划分。否则测试集里可能出现训练集中见过的相似抗体导致评估结果虚高。我通常按CDR-H3长度分组每组内再按序列相似度划分确保测试集和训练集之间的序列相似度低于70%。数据量方面目前公开的抗体结构大概几千个去掉冗余和低质量后能用的可能只有一两千。这个量级对于训练GNN来说偏少建议配合数据增强比如随机旋转、加噪声、或者用同源建模生成额外的结构。3.2 模型实现关键模块的代码骨架下面是一个简化的模型骨架用PyTorch和PyTorch Geometric实现。注意这只是核心逻辑实际代码还需要处理很多边界情况。import torch import torch.nn as nn from torch_geometric.nn import MessagePassing class AntibodyGNNLayer(MessagePassing): def __init__(self, node_dim, edge_dim, hidden_dim): super().__init__(aggrmean) self.node_mlp nn.Sequential( nn.Linear(node_dim hidden_dim, hidden_dim), nn.ReLU(), nn.Linear(hidden_dim, hidden_dim) ) self.edge_mlp nn.Sequential( nn.Linear(2 * node_dim edge_dim, hidden_dim), nn.ReLU(), nn.Linear(hidden_dim, hidden_dim) ) self.update_mlp nn.Sequential( nn.Linear(node_dim hidden_dim, hidden_dim), nn.ReLU(), nn.Linear(hidden_dim, node_dim) ) def forward(self, x, edge_index, edge_attr): return self.propagate(edge_index, xx, edge_attredge_attr) def message(self, x_i, x_j, edge_attr): return self.edge_mlp(torch.cat([x_i, x_j, edge_attr], dim-1)) def update(self, aggr_out, x): return self.update_mlp(torch.cat([x, aggr_out], dim-1))这个骨架里message函数负责计算每条边的消息update函数负责更新节点特征。实际使用时需要堆叠多层这样的GNN层并且在每层之后加入序列预测头和坐标预测头。序列预测头通常是一个线性层输出20维的logits。坐标预测头输出3维的偏移量。两个头共享GNN提取的特征但各自有独立的参数。3.3 训练技巧学习率、批大小与梯度裁剪训练GNN做抗体设计有几个超参数特别敏感。学习率方面我建议用1e-4到5e-4之间的值配合余弦退火调度。太高的学习率会导致训练不稳定损失震荡太低则收敛太慢。如果用了预训练初始化学习率可以再低一点1e-5到5e-5。批大小受限于显存。抗体结构通常有200-300个残基每个残基的特征维度在128-256之间。在单张24GB显存的卡上批大小大概能到8-16。如果显存不够可以用梯度累积来模拟更大的批大小。梯度裁剪是必须的。GNN在训练初期很容易出现梯度爆炸尤其是当图中有高度连接的节点时。我通常设裁剪阈值为1.0如果训练不稳定可以降到0.5。还有一个容易被忽略的点Dropout的位置。在GNN层之间加Dropoutp0.1-0.3有助于防止过拟合但不要在消息传递内部加否则会破坏图结构的完整性。3.4 评估指标除了RMSD和序列恢复率还要看什么评估生成质量不能只看两个指标。我一般会从以下几个维度综合评估指标含义理想范围序列恢复率生成序列与真实序列的匹配度40% (CDR-H3)结构RMSD生成坐标与真实坐标的偏差2.0 埃几何合理性键长键角异常的比例1%二硫键准确率正确形成的二硫键比例90%疏水补丁面积表面疏水区域大小100 平方埃多样性生成样本之间的序列差异熵 2.0特别要强调的是多样性。如果模型只会生成少数几种安全的序列即使单项指标很好看实际价值也不大。我通常会在测试时生成100-200个样本然后计算它们之间的平均成对序列差异。如果差异太小说明模型陷入了模式崩溃需要调整损失函数或增加采样温度。4. 实操中踩过的坑与排查思路4.1 生成结构出现原子重叠的根因定位训练到中期的时候我遇到过生成的结构里出现原子重叠的问题——两个残基的坐标距离小于1埃物理上完全不可能。一开始以为是学习率太高调低之后问题依旧。排查过程是这样的先检查了损失函数的几何约束项发现权重确实设得偏低只有0.05。把权重提高到0.3之后重叠问题明显减少但没有完全消失。然后我检查了消息传递的邻居定义发现只用了距离阈值导致一些在序列上远离但在空间上接近的残基没有被正确建模。改成序列邻近空间邻近的并集后问题基本解决。这个坑的教训是几何约束的权重不能太低而且邻居定义要同时考虑序列和空间两个维度。另外训练过程中要定期可视化生成的结构不要只看损失曲线。4.2 序列多样性崩塌模式崩溃的识别与缓解另一个常见问题是模式崩溃。模型生成的序列高度相似多样性极低。我遇到过一种情况模型对所有CDR-H3都生成相同的序列模式只是长度略有不同。识别这个问题很简单生成100个样本计算两两之间的序列一致性。如果平均一致性超过80%基本可以确定是模式崩溃。缓解方法有几种。一是增加采样温度让序列采样时更随机。但温度太高会牺牲序列合理性需要权衡。二是在损失函数里加一个多样性正则项惩罚生成样本之间的过度相似。三是用多个模型做集成每个模型用不同的随机种子训练然后合并生成结果。我试过最有效的方法是第二种在损失里加一个基于核密度的多样性项。具体来说计算当前批次内生成序列的成对相似度矩阵然后惩罚高相似度的对。这个项的权重设在0.1左右比较合适太高会导致生成质量下降。4.3 训练不收敛从数据到超参数的逐层排查训练不收敛是另一个高频问题。表现是损失震荡或者持续上升。我的排查顺序是这样的先看数据。检查输入的特征有没有异常值比如坐标里有NaN或者序列编码越界。这一步能排除大部分低级错误。再看模型。检查GNN层的梯度是否正常回传。有时候因为图结构太稀疏某些节点的梯度会消失。可以在GNN层之间加残差连接来缓解。然后看超参数。学习率是不是太大批大小是不是太小梯度裁剪阈值是不是设得太高我一般会做一个小的超参数扫描学习率试[1e-5, 5e-5, 1e-4, 5e-4]批大小试[4, 8, 16]梯度裁剪试[0.5, 1.0, 2.0]。最后看损失函数。各项损失的权重是否平衡如果某一项损失远大于其他项模型会被它主导忽略其他目标。我通常会在训练初期打印各项损失的值确保它们在同一量级。4.4 显存不够用时的几种实用降级方案抗体结构图比较大显存很容易不够。我试过几种降级方案按推荐程度排序第一种是梯度累积。用小的批大小比如2或4累积几步之后再更新参数。这能模拟大批大小的效果但训练速度会慢一些。第二种是子图采样。每次训练时只采样图的一部分节点和边而不是用完整的图。这能大幅降低显存占用但会引入噪声需要更多的训练轮数来补偿。第三种是混合精度训练。用FP16代替FP32显存占用能减少将近一半。但要注意GNN的某些操作在FP16下可能不稳定需要用动态损失缩放。第四种是模型并行。把GNN的不同层放在不同的GPU上。这需要改代码复杂度较高只在前面几种方法都不够用时才考虑。我一般优先用梯度累积混合精度这两个组合起来能解决大部分显存问题。5. 生成结果的后处理与实验验证5.1 结构弛豫让生成的结构更符合物理规律GNN生成的结构虽然在拓扑上合理但局部几何细节往往不够精确。比如键长可能有0.1-0.2埃的偏差键角可能有几度的偏差。这些偏差在后续的分子动力学模拟或对接实验中会导致问题。结构弛豫是必要的后处理步骤。我通常用两种方法一种是基于力场的能量最小化比如用AMBER或CHARMM的力场跑几百步最速下降另一种是基于知识的统计势比如用Rosetta的弛豫协议。实测下来能量最小化能把几何异常的比例从5%降到1%以下RMSD通常只会增加0.1-0.3埃是可以接受的。但要注意弛豫不能过度否则会把GNN学到的合理构象也破坏掉。我一般设置最大迭代步数为500步能量收敛阈值设在0.01 kcal/mol。5.2 结合对接实验验证结合亲和力生成抗体的最终目的是要能结合靶标。所以除了结构合理性还需要评估结合亲和力。常用的方法是对接docking比如用AutoDock Vina或HADDOCK把生成的抗体和靶标做刚性或半柔性对接。对接的流程大致是先从生成的抗体结构中提取CDR区域然后定义结合界面通常是CDR-H3和CDR-H1/H2再用对接工具搜索可能的结合构象最后用打分函数评估结合强度。这里有个坑对接的打分函数和GNN的训练目标不一致。GNN优化的是序列恢复率和结构RMSD而对接优化的是结合自由能。两者之间没有直接对应关系。所以GNN生成的高分样本对接打分不一定高。我的做法是生成大量样本比如500-1000个然后从中筛选出对接打分最好的几十个再做进一步分析。5.3 湿实验验证前的计算筛选策略在送去做湿实验之前还需要做一轮计算筛选排除明显不可行的候选。筛选维度包括聚集倾向用序列和结构特征预测抗体的聚集倾向。疏水补丁大、净电荷低的抗体容易聚集。免疫原性预测抗体序列中是否含有已知的免疫原性表位。表达产量预测抗体在常见表达系统如CHO细胞中的表达产量。稳定性预测抗体的热稳定性和pH稳定性。这些预测可以用现有的生物信息学工具完成比如用NetMHC预测免疫原性用Aggrescan预测聚集倾向。筛选之后通常能从几百个候选里留下十几个进入湿实验。我的经验是计算筛选能排除掉70%-80%的明显不可行候选大幅降低湿实验的成本。但筛选标准不能太严否则可能把真正有效的候选也排除掉。我一般会保留一些边缘候选用实验来验证。5.4 迭代精修轮数的选择太多还是太少都不行迭代精修轮数是一个需要仔细调的超参数。轮数太少结构没有充分优化质量不高轮数太多容易过拟合而且计算成本高。我做过一组对比实验在相同的初始化下分别跑10轮、20轮、30轮、50轮然后评估各项指标。结果是20-30轮时各项指标基本达到最优继续增加到50轮序列恢复率反而略有下降过拟合结构RMSD基本不变。所以我的建议是默认设30轮然后根据验证集的表现微调。如果验证集指标在20轮左右就饱和了可以降到20轮如果还在上升可以加到40轮。关键是要监控验证集指标而不是只看训练集。另外可以引入早停机制如果连续5轮验证集指标没有提升就停止迭代。这能节省不少计算时间。6. 这套方法目前的边界与我的个人体会迭代精修GNN做抗体序列-结构联合设计目前的能力边界还是比较清晰的。对于CDR-H3长度在8-15个残基之间的抗体生成质量比较好序列恢复率能到40%-50%结构RMSD能控制在2埃以内。但对于特别长的CDR-H3超过18个残基或者含有非标准构象的抗体生成质量明显下降。另一个限制是训练数据的量级。目前公开的抗体结构数据还是偏少模型容易过拟合。我试过用同源建模生成额外的训练数据有一定效果但提升有限。未来如果有更多高质量的结构数据这套方法的性能应该还能再上一个台阶。从实操角度我最大的体会是不要指望一次训练就能得到完美结果。抗体设计本身就是一个迭代的过程计算方法和实验验证需要紧密配合。GNN生成的结果更像是一个高质量的起点而不是最终答案。你需要用它来缩小搜索空间然后结合对接、分子动力学、以及湿实验逐步优化。最后分享一个小技巧在训练时可以先用一个小的、高质量的抗体子集做预训练然后再在完整数据集上微调。这样能显著加快收敛速度而且最终性能通常比直接在大数据集上训练要好。我试过用50个高分辨率抗体结构做预训练然后再用全部数据微调收敛轮数减少了约40%最终序列恢复率提高了3-5个百分点。