做 RosettaLigand 之前我每次都会在心里给“配体准备”单独留出半天。原因是它的报错模式太有迷惑性了多半不在准备这一步报而是等你跑完一轮 docking、看了三百个模型之后才在某个小角落里发现电荷不对、构象全是高能构象或者所有对接 pose 都叠在同一个地方。到了那一步再回头重新准备配体心态基本就崩了。这篇是“Rosetta 基础”系列的第二篇主题很窄——preparing ligand。不聊蛋白准备也不聊 docking 的 mover 怎么配只把“一个分子怎么变成 Rosetta 能认的 ligand”这件事讲透。你从这里能带走一套可以照抄的工作流从 SMILES 或 SDF 出发生成 3D 结构、构象库再通过molfile_to_params.py转成 params 和 PDB最后做一轮自检。刚接触 RosettaLigand、被 params 文件折磨过的人应该都能用得上。1. 配体准备为什么总是拦路虎理解 Rosetta 的配体表示方式1.1 一个分子在 Rosetta 里的两套身份很多新手会默认“配体就是 PDB 里的一段 HETATM”所以在prepare_receptor里检出了配体就觉得万事大吉。但在 Rosetta 的世界里配体远不止坐标这么简单。它需要一份 params 文件来描述原子类型、部分电荷、成键关系、内坐标甚至可旋转键和环的柔性。PDB 坐标只是身份卡params 才是它的完整简历。为什么 Rosetta 搞这么复杂因为打分函数要计算范德华、静电、溶剂化和氢键每个原子都必须有明确的原子类型和电荷采样时要调整可旋转键和环的构象程序必须知道哪些键能转、哪些原子是刚性的。没有这些信息Rosetta 没法对配体做 pack、做 minimize更别说采样了。这也是为什么直接拿一个普通 PDB 配体去跑复杂协议经常会报出一些让人摸不着头脑的错误——本质上就是这份“简历”没写清楚。1.2 你最终要拿到的三样东西跑完配体准备流程桌面上的关键产出有三个lig.params残基模板文件包含原子类型、电荷、成键、内坐标、可旋转键等所有参数。lig_0001.pdb配体的三维坐标文件。如果生成多个构象这个文件里会有多个MODEL ... ENDMDL块作为 Rosetta 的配体构象库。可选一个带上额外扭转信息的 SDF 文件某些 docking 流程里会用到。params 文件里有一行PDB_ROTAMERS lig_0001.pdb指向坐标文件。Rosetta 运行时会读取这行把配体的所有构象加载成 residue type 里的 rotamers。所以这两份文件必须放在同一目录下名字也得对得上否则后面报错会非常快。1.3 MOL2/SDF 到 params 的转换原理molfile_to_params.py做的事可以概括成三步读取输入分子的原子连接关系给每个原子分配 Rosetta 原子类型和部分电荷根据输入坐标计算键长、键角、二面角生成理想化的内坐标再把分子里的柔性自由度可旋转键、环等提取出来配合构象文件写成 params。本质上它是在把通用分子力场里的参数映射到 Rosetta 的残基模板语言。所以输入文件的质量会直接影响映射效果——原子类型标错、氢原子缺失、手性不对生成出来的 params 就是带病的。这个病不会在转换时报错而是在 docking 阶段以“很奇怪”的方式暴露出来。2. 从化学结构到 3D 坐标第一步不能省2.1 三种获取初始坐标的路线配体准备的起点是有一份靠谱的三维坐标。常见路线有三条各有取舍。路线一从实验结构中提取。如果你已经有 PDB 复合物直接把 HETATM 那一段抽出来当作初始坐标。好处是结合口袋里的真实构象就在眼前但缺陷是晶体结构常常有缺失原子、B 因子过高或者质子化状态根本不满足生理 pH。路线二从 SMILES 生成。这是最可控的方式尤其适合你只是想对一个已知化合物做 docking。用 RDKit 或者 OpenEye 的 Omega 都能从一维线性式构建合理的三维结构。下面这段是 RDKit 生成初始构象的标准操作from rdkit import Chem from rdkit.Chem import AllChem mol Chem.MolFromSmiles(CC(O)OC1CCCCC1C(O)O) # 阿司匹林 mol Chem.AddHs(mol) AllChem.EmbedMolecule(mol, AllChem.ETKDGv3()) AllChem.MMFFOptimizeMolecule(mol, maxIters500) w Chem.SDWriter(aspirin_3d.sdf) w.write(mol) w.close()路线三从 PubChem / ZINC 下载 SDF。胜在省事但下载下来的文件同样要检查手性、氢原子和质子化状态。我的建议是无论哪条路线拿到初始坐标后都做一次 3D 可视化和 valence 检查确认分子没有断键、没有奇怪的共价连接。2.2 质子化状态是精细活儿这是配体准备里最容易被忽视、却影响全局的一步。RDKit 的AddHs只是按化合价补全氢原子它不会考虑 pH。比如含羧基的分子在中性 pH 下通常以去质子化的羧酸根形式存在含氨基的分子则往往是质子化的铵离子。如果你把阿司匹林的羧基保持成中性酸氢键网络和静电分布都会明显偏掉对接结果自然不可信。一个实际的检查方法是以 pH 7.4 为默认值手动确认分子里的酸性基团和碱性基团状态。羧酸、磺酸、四唑这类通常去质子化脂肪胺、芳香胺、脒类通常质子化。如果口袋里的残基环境比较特殊还要结合氢键互补来判断。这一步没有现成命令能全自动解决OpenEye 的 QuacPac/FixpKa 可以给建议但最终拍板的人还得是你自己。3. 构象库生成docking 精度和性能的平衡点3.1 为什么 Rosetta 不像 MD 那样连续采样分子动力学是在连续空间里让键长、键角、二面角一直变化而 RosettaLigand 的做法是离散的先预生成一批构象组成配体位姿库docking 时在这个库里挑、在某些可旋转键上做扰动。所以构象库的质量直接决定了采样空间的天花板。构象太少可能漏掉真正能落在口袋里的活性构象构象太多pack 阶段的组合爆炸会让计算时间成倍上涨。这就是为什么“生成 200 个构象”不等于“生成 1000 个更好”。我试过的经验是常规小分子 50–200 个构象足够配合 0.5 Å 的 RMSD 聚类阈值能把冗余构象压掉又不至于错过主要自由度。3.2 我常用的生成命令和参数如果你用的是 RDKit可以基于 SMILES 生成多样化的构象集合from rdkit import Chem from rdkit.Chem import AllChem mol Chem.MolFromSmiles(CC(O)OC1CCCCC1C(O)O) mol Chem.AddHs(mol) params AllChem.ETKDGv3() params.randomSeed 0xf00d params.pruneRmsThresh 0.5 cids AllChem.EmbedMultipleConfs(mol, numConfs100, paramsparams) AllChem.MMFFOptimizeMoleculeConfs(mol, maxIters500) w Chem.SDWriter(aspirin_confs.sdf) for conf in mol.GetConformers(): w.write(mol, confIdconf.GetId()) w.close()如果你有 OpenEye 授权Omega 更方便omega2 -in aspirin.smi -out aspirin_confs.sdf \ -maxconf 100 -rms 0.5 -ewindow 10这两个命令背后都在做一件事在能量窗口范围内找多样化的低能构象。-ewindow 10表示只保留相对最低能量 10 kcal/mol 以内的构象-rms 0.5则是去掉空间叠合度太高的冗余结构。这两个参数值得你根据分子大小调一调环多、柔性大的分子要放宽能量窗口扁平芳香分子则可以收紧。4. molfile_to_params.py 全流程实战拆解4.1 一条完整命令构象文件准备好之后进入核心转换环节。以阿司匹林为例完整命令如下python3 $ROSETTA/main/source/scripts/python/public/molfile_to_params.py \ -n LIG \ -p lig \ --conformers-in-one-file \ aspirin_confs.sdf解释一下参数-n LIG是给配体残基起的三个字母名字这个名字必须和复合物 PDB 里的 HETATM 残基名一致-p lig是输出文件名的前缀会生成lig.params和lig_0001.pdb--conformers-in-one-file表示把所有构象写进同一个 PDB 文件每个构象作为一个MODEL。如果不用这个参数脚本会按_0001、_0002生成多个 PDB 文件后面指定起来麻烦很多。这里有个容易踩的细节老版本 Rosetta 的脚本可能需要python2新版本基本都在 Python 3 环境跑。报错信息里如果出现SyntaxError先别怀疑自己的命令先确认解释器版本对不对。4.2 输出文件里到底有什么打开lig.params你会看到类似下面这样的结构NAME LIG IO_STRING LIG Z AA UNK TYPE LIGAND ATOM C1 CH3 X 0.153 ATOM C2 CA X 0.021 ... BOND C1 C2 BOND C2 C3 ... IC C1 C2 C3 C4 1.523 120.0 180.0 120.0 1.420 NBR_ATOM C3 NBR_RADIUS 5.2 PDB_ROTAMERS lig_0001.pdbATOM行里依次是原子名、Rosetta 原子类型、电荷和内部编号。BOND和IC行定义了成键关系和内坐标。PDB_ROTAMERS指向构象文件。如果你单独把lig_0001.pdb拷走、把 params 留在原目录Rosetta 在加载时会找不到坐标文件这个坑我见过不止一次。到这一步你其实已经可以把这个 ligand 放进 RosettaLigand 了运行任何相关协议时加上-extra_res_fa lig.params即可。如果在复杂体系里也可以用-in:file:extra_res_fa一次性传多个 params。做这一步的意义是让 Rosetta 把配体当“真正的残基”读进去而不是当成一段没人管的杂原子。4.3 能量最小化要不要做molfile_to_params.py有一个--minimize选项会在转换前对输入构象做一次局部优化适合那些从 2D 绘制的分子、初始坐标本身比较粗糙的情况。如果你输入的是已经经过 MMFF 或 Omega 优化的构象--minimize不是必须的。我的习惯是转出来的 params 先用后继的 scoring 测试跑一遍分数异常才回头重转。节省的时间远大于省下的麻烦。5. 电荷、原子类型、手性最容易被忽略的三个隐形地雷5.1 电荷的准确性直接影响打分Rosetta 的ref2015打分函数里静电能和溶剂化能都跟原子电荷强相关尤其是带电荷的配体。molfile_to_params.py默认会优先用 OpenEye 的 AM1-BCC 电荷如果环境里没有 OpenEye 会退回 Gasteiger。AM1-BCC 对极性和氢键供受体的描述更细Gasteiger 则是老牌快速方法胜在开箱即用。不管用哪种我都建议你做一个反直觉但极其重要的自检手动把 params 文件里所有ATOM行的电荷加起来它应该等于这个分子在目标质子化状态下的形式电荷。羧酸根是 -1铵离子是 1中性分子是 0。这一步如果对不上后面所有静电相关的结果都不值得看尽早排查。5.2 原子类型与命名原子名是 Rosetta 里一个很“较真”的地方。ATOM名必须与 PDB 文件里的原子名一致任何一个错位都会导致内坐标计算错误。转换脚本通常会自己分配原子名但如果输入 SDF 里的原子名重复、过长或包含非法字符就可能在后面步骤里出问题。遇到这种情况可以先在初始 SDF/MOL2 里把所有原子名改成唯一的、不超过三四个字符的名字再重新转换。原子类型识别失败也比较常见尤其是含卤素、金属、膦酸等不太常规的基团。Rosetta 的原子类型体系是从化学环境推出来的输入文件的格式越标准识别越稳定。我个人的建议是优先用带氢原子的 SDFMOL2 次之最怕的是那种连氢都没有、靠脚本猜价态的“光杆”结构。5.3 手性不对对接方向直接反掉这个坑比前两个更隐蔽。从 SMILES 生成 3D 结构时如果 SMILES 里没有明确的或手性标记生成器可能给你任意一个对映异构体。Rosetta 不会检查“这个是不是我想要的构型”它只负责按你给的坐标算能量。你把 R 构型当成 S 构型拿去对接最后出来的优选 pose 很可能就是错误的镜像结合模式。所以拿到初始 3D 结构后务必要用你知道的手性中心逐个确认。有实验结构就优先参考实验结构没有实验结构至少要在 PyMOL 里把分子转一圈看看关键手性中心的构型是否符合预期。这一步花不了几分钟却能避免整轮 docking 白跑。6. 常见报错对照与一次完整的自检6.1 高频报错与含义我把这些年见过的高频问题整理成了一张表基本都是配体准备阶段可以现场排查的报错或现象常见原因处理思路Unable to find params for LIG没加-extra_res_fa或路径不对检查命令行里是否带了 params 文件路径用绝对路径Unrecognized residuePDB 里的 HETATM 残基名与-n参数不一致统一三个字母的残基名PDB_ROTAMERS file not foundparams 和 PDB 没放在同目录把lig_0001.pdb和lig.params放在一起再跑某个原子没有合法原子类型输入文件原子类型、杂化状态有问题换成 SDF 并显式带上氢原子检查是否有非标准元素对接分数高得离谱电荷异常或质子化状态错误检查总电荷、pH 状态必要时用--minimize重转构象库只读到一个构象输入 SDF 只有一个构象或没加--conformers-in-one-file确认 SDF 里构象数量并用多构象 SDF 重跑6.2 上考场前我必做的四步自检配体准备完别急着跑完整 protocol。我的习惯是先做四步快检。第一步跑一次 scoring 冒烟测试$ROSETTA/main/source/bin/score_jd2.default.linuxgccrelease \ -s lig_0001.pdb \ -extra_res_fa lig.params \ -out:nooutput没有报错且分数在合理范围不是 NaN、不是几百上千说明 params 至少能被 Rosetta 正常加载。第二步验证总电荷。把 params 里所有ATOM行的电荷手动加和对不上就立刻回头检查质子化状态。第三步用 PyMOL 可视化lig_0001.pdb。重点看键长键角是否合理、有没有原子重叠、手性中心是否正常。这一步对新手尤其重要很多问题靠眼睛就能发现。第四步检查可旋转键数量。在 params 文件里数一下你关注的柔性二面角和分子结构做对比。常见小分子的可旋转键数量一般是个位数如果多出很多说明脚本可能把某些刚性键也当成可转键了这时候需要回到输入文件的原子类型去排查。四步都过了这个 ligand 才算真正能用。我在实际操作中见过太多项目卡在 docking 打完才发现配体没有准备干净重新跑一轮的时间成本远远高于一开始认认真真做这十几分钟自检。Rosetta 不会替你判断配体准没准备好它能做的只是把你喂进去的东西当成事实来尊重——所以你喂进去的每一步都得自己心里有数。