中药靶点筛选计算实操（二） | 分子对接：针对明确靶点，从中药库中筛选高置信先导并预测结合模式

当反向对接或文献/组学数据已经给出了一个或多个明确的候选靶点，我的研究目标就从“找靶点”切换为“从中药成分库中筛选能与该靶点高效结合的分子，并解析其具体结合模式”。此时使用正向分子对接。

核心任务是：以该靶点蛋白的三维结构为受体，以中药化合物库（成千上万个）为配体，通过打分函数排序，选出打分靠前的候选物，并观察其与关键残基的相互作用。

从PDB下载靶点结构（优先选择分辨率<2.5 Å、有共结晶配体的结构）。随后执行以下步骤，全部在Schrödinger Suite或MGLTools中完成：

除少数被文献证实的“结构水”（参与氢键桥接）外，其余水分子全部删除。我使用Protein Preparation Wizard中的“Delete waters beyond 5 Å from ligand”选项，再手动保留与配体形成双齿氢键的水。

将pH设为7.0~7.4，使用PROPKA计算可滴定残基的pKa。特别注意组氨酸（His）——根据局部环境判断是HID（δ-质子化）、HIE（ε-质子化）还是HSP（双质子化）。这一步骤若出错，后续氢键网络将完全失真。

仅对添加的氢原子进行约束性最小化（OPLS3e力场，RMSD收敛阈值0.3 Å），重原子坐标固定不动，避免骨架偏离原始晶体结构。

我的中药化合物库包含TCMSP中所有口服生物利用度≥30%且类药性≥0.18的分子（约4000个）。每个分子经历：

❶ 从SDF或SMILES生成三维构象（使用LigPrep或Open Babel，生成最多32个立体异构体和互变异构体）。

❷ 分配Gasteiger电荷，指定可旋转键。

❸ 导出为PDBQT（用于AutoDock Vina）或MAE（用于Schrödinger Glide）。

各维度从中心向外扩展至覆盖口袋的全部残基，再额外增加5 Å。例如，若口袋半径为10 Å，则盒子边长设为30 Å（即±15 Å）。

设为0.375 Å（默认），对于口袋狭小的靶点可降至0.25 Å，但会显著增加计算时间。

对于初筛，exhaustiveness=8已足够；对于最终确认的少数候选物，我提升至20。

我偶尔需要筛选包含天然产物衍生物的大型库（如ZINC的天然产物子集，约5万分子）。此时直接用Vina耗时过长，采用串联策略：

使用基于分子指纹的深度全连接网络（DFCNN）模型进行快速分类。我使用开源DeepBindBC或自己训练的模型，输入为ECFP4指纹，输出为“结合概率”二分类。设置阈值，保留概率>0.8的分子（通常滤掉80%）。

仅对保留的约20%分子进行Vina精确对接。这使总计算时间从数天缩短至数小时。

检查是否配体穿透了蛋白表面（即原子间碰撞距离<1.5 Å）。我使用PyMOL的“clash”检测功能，若有严重碰撞，舍弃该构象。

以共晶配体的结合姿态为标准，计算RMSD。若RMSD>2.5 Å，即使打分好也存疑，往往是力场参数缺陷或构象采样不足导致。

在PyMOL中显示氢键，确保供体-受体角度>120°。角度过小则视为弱相互作用，不计入有效氢键。

正向分子对接完成了“哪个分子以什么姿势结合靶点”的预测。它的输出是静态的 —— 没有溶剂，没有热运动，没有时间演化。

在我的经验中，对接得分最高的分子在后续实验中失败的比例依然很高，因为许多高得分构象在真实水环境中会迅速失稳。因此，我只将对接结果作为“候选列表”的最终筛选器，通常保留排名前5~10的化合物。

但若要真正判断这些候选物能否在纳秒级动态环境中维持结合，必须进入分子动力学模拟。

下一篇，我将详细讲述如何运行100 ns MD模拟，以及如何用gmx_MMPBSA计算结合自由能，从而剔除动态假阳性。敬请期待后续实操内容。

点击文末阅读原文，

获取更多中药作用机制研究实验方案~

同时，如果您希望将上述计算流程快速整合到自己的课题中，或需要从靶点发现到分子验证的一站式支持，欢迎了解以下专业服务——

参考资料