干过分子动力学模拟的人都知道蛋白结构准备这步看着不起眼实际决定了一个体系能不能顺利跑起来。我在刚入门Amber那阵拿着一个从PDB下载的晶体结构改改残基名就直接load进去结果tleap一路报错问了师兄才知道光是“读入结构”这一步就有一堆门道序列对不对齐、质子化状态有没有确定、缺失残基要不要补、二硫键该不该显式声明……任何一个环节马虎后面三五十纳秒的模拟都可能是白跑。这篇文章就把我整理出来的Amber蛋白结构准备流程和关键注意事项完整讲一遍适合刚接触Amber、准备拿蛋白体系做模拟的师弟师妹也适合那些已经能跑通最小化、但经常在结构上翻车的朋友。1. 结构准备没做好后面全是无用功1.1 我踩过的结构准备坑刚做Amber那会儿我拿到一个膜蛋白的冷冻电镜结构心想电镜分辨率也够于是我直接tleap loadpdb结果报“Unknown residue: MSE”。我以为是电镜结构的问题后来才明白MSE是硒代甲硫氨酸的PDB命名Amber标准残基库认不了要先转成MET或者用特殊参数。这不是什么深奥的原理问题纯粹是结构准备阶段该处理的事。还有一次更隐蔽体系净电荷算出来是6我顺手加了6个Cl-结果最小化的时候体系一直在膨胀后来检查发现蛋白里还有两个钙离子我没删净电荷其实是14离子浓度完全不对。这些事故其实都可以在结构准备阶段通过简单的检查避免但如果你把结构准备当成“导入即用”的步骤就很容易踩进去。1.2 结构准备到底在准备什么从生物学的角度讲PDB库里的蛋白是实验测得的静态构象而分子动力学模拟要的是这个蛋白在特定pH、盐浓度、温度下的动态行为。结构准备的核心工作就是把实验结构转换为一个“在给定模拟条件下热力学合理”的初始状态。具体来说包括这五件事处理实验结构中的“杂质”水、配体、离子、去污剂哪些该留哪些该删补充缺失的原子和残基保证肽链连续、侧链完整根据模拟环境的pH设置正确的质子化状态明确二硫键、金属配位等共价或配位相互作用用力场能识别的残基名和格式重建拓扑输出后续计算需要的参数文件。通过这些操作最终要得到两个核心文件prmtop拓扑文件和inpcrd坐标文件。prmtop包含原子类型、电荷、键连关系、力场参数、水盒子尺寸等所有分子力学参数inpcrd包含每个原子的初始三维坐标。后续能量最小化、升温、平衡、生产模拟都以这两个文件为起点。这两个文件一旦有错后面所有阶段的输出都不可信。1.3 为什么要把“一般流程”单独拿出来讲很多教程是从“加载PDB、加溶剂、加离子、跑起来”四步走开始的但实际项目中结构来源千奇百怪有的从X射线晶体衍射来有的从冷冻电镜来有的是同源建模拼出来的有的来自分子对接的复合物。不同来源的结构表面“长得像蛋白”内部却带着各自的遗留物。这篇文章讲的就是一套能应对这些不同来源的通用检查流程我把它拆成了五件事拿结构、清结构、定质子化、建拓扑、验拓扑。每一步都不难但顺序和细节决定了成败。2. 从PDB文件到“干净”的蛋白结构2.1 把PDB文件当成一本“账本”来看PDB文件看起来是一堆行文本但它是一个结构化的账本。拿到结构第一件事不是急着loadpdb而是先把文件“读一遍”。重点看这几个区块HEADER分子类型和来源能判断这是晶体结构、核磁结构还是电镜结构SEQRES完整的氨基酸序列包括晶体里看不到的残基ATOM蛋白质原子实际坐标这些是有实验电子密度/密度图支撑的部分HETATM非标准残基、配体、水等SSBOND二硫键记录REMARK 465缺失残基的说明晶体结构里经常有残基标了“缺失”或“没有电子密度”。我实际中最常做的第一步是用PyMOL或VMD打开结构显示成连线或卡通模式看序列是否完整、链编号是否正确。这不费时间但能提前发现很多脚本层看不到的问题。比如有的结构里同一个氨基酸出现两个交替构象双占位PDB里的ALT A/B标记PyMOL里看就是同一位置有两套原子如果不处理tleap就会因为同一个残基出现重复原子而报错。另外分辨率这个指标也必须看X射线晶体结构分辨率在2.5 Å以下比较可靠冷冻电镜看整体分辨率和局部密度3 Å左右是合格线同源建模结构看序列相似度和模板覆盖度。分辨率太差的结构后面模拟出来的局部构象可信度也有限应该在结构准备阶段就作出判断。2.2 缺失残基与不完整侧链补还是不补晶体结构里loop区或膜蛋白的膜外区域经常缺失残基PDB文件里REMARK 465会列出所有缺失的残基和原子。处理方式取决于缺失的大小缺失单侧链或少量的几个骨架原子通常直接用tleap补全Amber在构建拓扑时会自动加缺失原子这是最高效省事的方案缺失一段loop通常多于5个残基建议用同源建模工具如Modeller或Swiss-Model补全。不补而直接模拟会导致蛋白链断成两截拓扑错误。而且真正影响功能的关键残基可能就在缺失区域内比如酶的底物结合loop跳过它模拟意义就打了折扣。值得提醒的是补完的loop没有实验坐标约束它的构象是建模工具推出来的不代表真实状态。所以在模拟流程里建议对补环区域先做约束最小化或短时间约束平衡再放开全局模拟。否则一开始的全局最小化会让这段loop以相当人工的方式塌缩到蛋白表面之后的模拟结果很难说靠谱。2.3 用pdb4amber把格式标准化AmberTools自带pdb4amber脚本是结构准备的第一步标准化操作。它能把PDB里的命名转换成Amber能识别的形式还能去掉水、识别链、加氢。典型用法pdb4amber -i input.pdb -o protein_clean.pdb --reduce这个命令会生成几个文件protein_clean.pdb处理后的主文件protein_clean.pdb_sslink识别到的二硫键protein_clean.pdb_renum.txt残基编号映射。我现在的习惯是先不加--reduce跑一遍看看输出pdb4amber -i input.pdb -o protein_clean.pdb不加氢的版本只是清理之后在tleap里由力场自动加氢这对一般蛋白更省心。--reduce会调用reduce程序加氢对学生物出身的用户多了一层概念负担而且在有配体或非标准残基时reduce加出来的氢经常让人看不懂。如果非用--reduce千万记得把输出文件打开肉眼扫一遍特别是半胱氨酸巯基上的氢和二硫键附近别盲信自动处理。pdb4amber还支持-d参数直接删除水-y参数保留整个残基而不是把缺失原子标成“未知”。我经常组合成pdb4amber -i input.pdb -o protein_clean_nohyd.pdb -d删水这步放到这里做最顺手后面tleap里就少一步删水操作。2.4 水、配体、离子留还是不留从晶体结构里去掉水是有学问的。绝大多数结晶水在模拟中会与体相水交换删掉没影响。但位于活性位点、参与底物配位或稳定关键残基的保守水分子有时对功能模拟很重要。我的原则是先查文献如果文献强调某个水分子在催化或结合中的作用就保留它其余水全部删除。至于配体分三种情况课题不研究配体结合直接删除但要注意配体存在时引起的蛋白构象变化删掉后蛋白可能因为“空出来一个大空腔”而在模拟初期发生局部塌陷这属于正常现象不必紧张做蛋白-小分子复合物模拟配体要单独参数化。配体坐标从PDB提取成mol2或pdb格式在antechamber里用GAFF力场生成frcmod和lib文件这一步要走独立的流程配体是共价结合的抑制剂、底物类似物需要在力场参数里定义共价连接处理更复杂通常要用到amber的链接原子或者特殊的构建方式。离子方面钙离子、锌离子等金属离子在一些蛋白中起结构稳定作用需要保留。但保留离子意味着你得确认它的配位环境tleap不会自动为金属离子建立配位键。对结构性的钙离子一般的做法是在tleap里用它所在区域的水和蛋白残基建立非键模型或者直接用一些专门的金属中心参数。如果只是模拟缓冲液中的游离钠氯离子则完全交给addions处理。3. 质子化状态pH不是小事3.1 为什么必须显式加氢PDB实验结构几乎不提供氢原子坐标但在分子力学计算中氢原子的存在直接决定了残基的带电状态和氢键网络也就决定了静电相互作用。举个例子组氨酸咪唑环的Nδ和Nε哪个带氢直接决定His能不能作为质子供体或受体这在酶催化机制模拟里是成败关键的。Amber默认的tleap加氢逻辑会给每个残基补全在其标准质子化状态下的氢在pH7附近Asp和Glu去质子化带负电Lys和Arg质子化带正电His默认成中性但形式需要自己选择。如果你不额外处理跑的就是“标准pH 7.0”的蛋白这在多数情况下正确但有些残基出现了远离7的pKa就需要手动干预。3.2 组氨酸、半胱氨酸等特殊残基的质子化判断组氨酸在Amber残基库里有三种命名这个必须背下来HIDNδ带氢中性HIENε带氢中性HIP两个氮都带氢带1正电荷。对很多酶来说His到底是HID还是HIE可以通过氢键网络大致判断咪唑上的N如果能氢键给受体一般带氢如果作为氢键受体则不带。更严格的做法是算pKa用H或PropKa工具。半胱氨酸也有两种形式CYS硫醇带巯基氢中性CYX硫原子参与二硫键不能加氢中性。这里有一个特别容易翻车的点tleap不会自动识别二硫键。如果PDB里有SSBOND记录你需要显式地用bond命令把两个Cys的SG原子连起来同时把残基名改成CYX。如果不改两个半胱氨酸会各自保留巯基氢形成错误的硫醇状态如果只改名不建键则拓扑里会缺失关键的共价连接。3.3 实际判断工具H和PropKa现在用得比较多的是H服务器和PropKaH在线服务器输入PDB和pH、盐浓度输出质子化状态和pKaPropKa本地命令行速度快适用于批量处理。以H为例提交后它会给出每个可滴定残基的pKa和建议质子化状态按建议把PDB里的HIS改成HID/HIE/HIP把CYS改成CYS或CYX即可。我有一条个人经验H的输出通常偏保守如果某个残基的pKa卡在6.9这种离pH7很近的位置我会把两种状态都建出来分别跑几十ns看哪个更稳定而不是拍脑袋选一个。还有一点很多PDB文件里的残基是“旧式命名”比如组氨酸叫HSD或HSE而不是HID/HIE天冬酰胺和亮氨酸的原子命名也和老版本不完全一致。pdb4amber能自动把大部分旧命名转成新命名但最好还是对照转换后的PDB确认一遍。常有人问我为什么Amber教程里组氨酸要写HID自己的PDB里明明是HIS——因为PDB用HIS泛指组氨酸Amber则区分质子化微观态这是两套命名体系不是写错了。4. 用tleap构建拓扑和坐标脚本逐行拆解4.1 选力场ff14SB和ff19SB怎么挑对于普通可溶性蛋白ff14SB是当前最常用的选择它改进了侧链扭转角和骨架二面角的描述平衡性好文献支撑也多。ff19SB是对ff14SB的进一步更新对骨架的二面角参数做了系统性再拟合在某些体系上表现更好尤其对抗原-抗体这类对骨架构象敏感的体系。但我的建议是新手先从ff14SB起步。原因有三个第一ff14SB的参数验证时间长可迁移性强第二它配套参数齐全修饰氨基酸、配体、脂质等周边参数的兼容性好第三网上踩坑案例多出问题好搜索。如果做膜蛋白可选择ff14SB配合脂质分子力场做核酸则用OL15或bsc1做蛋白质-糖复合物则要确认糖部分的力场配套情况不是简单改一个leaprc就能解决的。tleap通过加载不同的leaprc文件进入不同力场环境source leaprc.protein.ff14SB source leaprc.water.tip3p这里还有个细节如果体系和金属离子相关可能需要加载leaprc.ff14SBplus或专门的金属参数文件如果体系里有特殊修饰还要加载对应的frcmod文件。总之力场是结构准备里“选参数”的一层所有选择的依据都应该记录到工作日志里。4.2 一个完整可用的tleap输入脚本下面是我常用的脚本适用于一个不含配体的可溶性蛋白# tleap.in source leaprc.protein.ff14SB source leaprc.water.tip3p source leaprc.gaff mol loadpdb protein_clean.pdb # 二硫键示例残基23和残基87 bond mol.23.SG mol.87.SG # 加溶剂盒子蛋白边缘距盒子边界10 Å solvateoct mol TIP3PBOX 10.0 # 加NaCl到0.15 M同时中和净电荷 addions mol Na 0.15 addions mol Cl- 0.15 # 保存 saveamberparm mol prmtop inpcrd savepdb mol check.pdb quit逐个解释一下关键语句source命令加载参数文件顺序不能乱先加载蛋白力场再加载水模型最后加载GAFF用于后面可能要处理的配体或特殊残基loadpdb读入标准化后的结构这一步tleap会自动补全缺失原子同时会输出一些信息比如“Added missing atoms”或“WARNING”需要仔细看bond mol.23.SG mol.87.SG建立二硫键这里mol.23表示蛋白的第23个残基SG是硫原子名solvateoct mol TIP3PBOX 10.0用TIP3P水模型加截角八面体盒子10.0是蛋白表面到盒壁的最短距离单位是Å。八面体盒子比正方体体积效率高做蛋白旋转扩散模拟时是标配addions mol Na 0.15这句同时做了两件事先中和体系净电荷再加到0.15 M离子浓度。所以脚本里先加Na再加Cl-顺序上一般不会错saveamberparm输出prmtop和inpcrd这是核心产出savepdb mol check.pdb额外输出一个检查用PDB方便在PyMOL里验证有没有原子重叠、离子位置是否合适。4.3 关于水盒子类型的选择tleap里有几种常见的溶剂添加方式solvatebox正方体盒子简单直接solvateoct截角八面体原子数更少边界效应更小我几乎所有可溶性蛋白都用它solvatedock或solvent cap球形溶剂帽只用于非周期边界条件的小体系、或者做局部采样的特殊情况。表面蛋白、球蛋白用solvateoct就好如果做薄膜体系或需要各向异性才考虑solvatebox。盒子大小方面经验值是10 Å为下限缓冲液太薄会引入周期性镜像相互作用太厚则白白增加原子数、拖慢计算。常用的是10~12 Å如果你做长链核酸或柔性较大的体系可以放宽到12~15 Å。看完prmtop里的BOX_SIZE能及时发现盒子设置是否合理。4.4 多个链、端基封端和修饰残基的处理如果蛋白是多聚体loadpdb一次载入即可但要注意残基编号冲突。tleap不关心CHAINS它按残基编号区分原子如果两个链的残基编号都是从1开始则会因为残基编号重复而报错。解决方法是先用pdb4amber的renum功能对第二个链重新编号或者在载入前用PyMOL脚本给不同链的残基重新编号。端基封端也常常被忽略。如果模拟的是一个截短的肽段或结构域的端基裸的带正电的N端和带负电的C端会使末端残基的静电行为不真实。常见做法是加ACE乙酰基封N端、NMEN-甲基酰胺封C端。在PDB里把N端残基前加一个ACE残基或把C端残基后加NMEpdb4amber能识别这些修饰名tleap也能直接载入。磷酸化、糖基化、甲基化等翻译后修饰在Amber里需要专门的参数。tleap能够识别部分磷酸化残基如TPO、SEP、PTR但前提是加载对应的frcmod参数文件。如果修饰不含Amber官方参数就要用类似GAFF参数化的小分子流程单独处理这已经是另一节课的内容了。简单提醒一句遇到修饰残基先查Amber的残基库列表再决定是直接加载还是走小分子参数化流程。5. 提交最小化前的结构与拓扑验证5.1 检查净电荷和整体合理性tleap跑完先别急着提交任务。用下面命令快速查看prmtop信息ambpdb -p prmtop -c inpcrd final.pdb或者在tleap脚本里打印残基数量、原子数量。用parmed也可以读取拓扑并查看电荷信息parmed prmtop # 在parmed里 # charge :* # 输出后把电荷求和净电荷不是0的话模拟中周期性盒子无法正确处理静电项sander/pmemd会在运行时给出严重警告或直接报错。经验做法是让净电荷为0或小整数电荷配合反离子中和。顺便把check.pdb在PyMOL里打开检查三件事蛋白质链是否有异常断点有没有原子落在盒子外面加进的离子有没有落到蛋白内部。如果离子落到蛋白内部通常是因为蛋白有较大空腔或通道属于正常现象但如果离子落在疏水核心附近就要警惕了说明离子浓度或初始摆放可能有问题。这种情况下可以尝试减少离子浓度或者改用addions时指定离子放置区域用around和exclude等参数。5.2 我整理了一份高频报错速查表下面这些是我实际跑tleap时碰到过的高频报错列成表格供参考报错信息根因处理办法Unknown residuePDB里有力场不认识的残基名比如MSE、HSD或配体用pdb4amber转换残基名或者为配体单独准备frcmod文件does not have a type某个原子找不到对应力场原子类型检查是否加载了正确的leaprc和frcmod确认非标准残基已处理Created a new bondtleap根据距离自动建了新键通常是初始构象中原子间距异常可视化检查该区域的坐标往往来自错误的残基编号或缺失残基补键WARNING: The unperturbed charge of the unit is not zero体系净电荷不为0用addions重新中和或检查是否漏删了未处理的离子最小化时体系能量爆炸初始结构中有原子重叠或坏接触检查是否有重叠原子必要时先用最陡下降法小步优化或给蛋白加位置约束先优化水这里特别说一下“Created a new bond”这条。tleap会根据原子间的距离自动判断是否成键。如果两个本不该成键的原子距离过近它会警告并自动建键。这种情况十有八九是坐标文件里有重叠原子或者残基编号错位导致某个残基的骨架和另一个残基的骨架被算成了同一个残基。遇到这条警告不要忽略一定要找到生成新键的区域可视化后确认是哪种原因。5.3 用一次成功的能量最小化验证准备结果验证结构准备是否合格最快的方法是跑一个短小的能量最小化。最小化阶段如果能量能稳定下降并收敛说明拓扑和初始坐标基本合格。一个常用的sander/pmemd的min.in模板如下cntrl imin1, maxcyc5000, ncyc2500, cut10.0, ntb1, ntc1, ntf1, ntpr100, /解释一下关键参数imin1表示做最小化ncyc2500指定先做2500步最陡下降法Steepest Descent。最陡下降法处理坏接触的效率高但它收敛慢所以跑2500步后自动切换成共轭梯度法直到maxcyc5000步为止cut10.0是非键相互作用的截断距离单位ÅAmber中一般用10ntb1表示在周期性边界条件下做最小化ntc1和ntf1表示不约束键长这对最小化阶段是合适的。有些教程让用户用ntc2约束涉及氢的键反而会在初始构象有问题时掩盖坐标异常。运行命令pmemd.cuda -O -i min.in -p prmtop -c inpcrd -o min.out -r min.rst -ref inpcrd如果没有GPU把pmemd.cuda换成pmemd或sander.MPI即可。跑完后看min.out最后几行的ENERGY值重点关注GRMS均方根梯度和总能量。GRMS在0.1 kcal/mol·Å以下总能量稳定在一个平台区间就基本可以认为结构准备这一关过了。很多教程直接跳到平衡和生产我认为最小化这一小步无论如何不要省。它的计算代价极低却是最廉价的整体自检手段。一套拓扑如果有问题最小化阶段往往会以能量爆炸、原子飞出的方式暴露出来如果最小化能顺利收敛后面跑平衡和生产的安全感会强很多。6. 最后再分享几个实用细节这篇写到这里把蛋白结构准备的主流程遮了一遍。最后补充几个我实际用下来的细节都是踩过之后才长记性的。养成记录每次结构准备操作的习惯。哪个PDB条目、去了哪些水、补了哪些残基、电荷怎么中和、用的哪个力场和版本全都写进README或日志。等到模拟结果异常回查起来要省太多时间。我现在每个项目目录里都有一份prepare_structure.log一行行记操作和时间出问题排查效率高得多。pdb4amber的输出文件里有一些看名字容易忽略的附属文件比如_sslink和_renum.txt。前者记录识别到的二硫键后者是残基编号映射。强烈建议每次都打开看一眼特别是从多聚体、有替代构象的结构出发时。我遇到过二硫键在pdb4amber里识别对了但我自己在处理PDB时不小心删了相邻的残基导致SSBOND记录和实际坐标对不上tleap里建键就报错。拿不准的质子化状态多建几版体系各跑一个最小化加短平衡比较能量和结构波动再决定。这个成本比生产模拟低很多却能把不确定因素控制在源头。一个小技巧是看这几版平衡后His周边氢键网络是否和文献描述一致如果不一致大概率是质子化状态选错了。结构准备的很多“标准答案”要结合实验背景来理解你模拟的是一个酶活pH为4的蛋白却按pH 7做了标准质子化跑再久也解释不了实验数据。任何工具给出的质子化建议都要回到实验条件里验证一遍。实际上Amber模拟的成败并不取决于生产相位跑了多长、采样有多广多数问题在结构准备这关就已经埋下了。把结构准备当作一个独立的、需要反复检查的流程来对待而不是一个开机即用的导入步骤后面所有阶段的顺利程度都会明显不一样。我这里给出的是一套通用流程具体到每个蛋白、每个突变、每个配体体系还会有各自的细节要补充到时候再逐个展开聊。