gnina分子对接实战：CNN打分从安装到批量筛选

1. 分子对接的现状与gnina的破局点做药物设计或者计算化学的朋友对分子对接这四个字肯定不陌生。不管是做虚拟筛选还是研究蛋白-配体相互作用模式分子对接都是绕不开的一步。过去十几年大家用得最多的无非就是那几个商业软件比如薛定谔的Glide、MOE的Dock、或者开源阵营里的AutoDock Vina。这些工具各有各的脾气商业软件界面友好、流程成熟但价格不菲开源工具免费但打分函数的老化问题一直被人诟病。gnina的出现算是给这个领域投下了一颗不大不小的石子。它的核心卖点很直接用卷积神经网络CNN来做分子对接的打分。传统打分函数不管是基于力场的还是基于经验的本质上都是手工设计特征的线性组合或者简单非线性组合。而CNN不一样它能自动从原子的空间分布中学习高维特征理论上能捕捉到传统打分函数漏掉的复杂相互作用模式。我第一次接触gnina是在一个虚拟筛选项目里当时用Vina跑了一遍命中率不太理想假阳性偏高。后来看到gnina的论文说是在多个基准测试上比Vina和Vinardo都好就决定试一试。实测下来CNN打分版确实在富集因子EF上有明显提升尤其是早期富集。这篇文章就把我从安装到跑通第一个案例的完整过程拆开来讲包括踩过的坑和后来总结出来的技巧。gnina本身是基于AutoDock Vina的代码框架改的所以如果你之前用过Vina上手会很快。但它又不仅仅是换了个打分函数它把CNN打分和传统的Vina打分做了融合还支持多种对接模式。最关键的是它完全开源代码在GitHub上可以自由获取。对于预算有限但又想做高质量对接的课题组来说这几乎是最优解。不过gnina的安装和配置比Vina要复杂一些因为它依赖libmolgrid这个库来做分子网格化还需要OpenBabel来处理分子格式转换。这些依赖关系如果没理顺编译的时候会报一堆错。我见过不少人在安装这一步就放弃了转回去用Vina。其实只要把依赖装对gnina的编译并不难。下面我就按实际操作的顺序把整个流程讲清楚。2. 环境准备与依赖安装的实操细节2.1 系统选择与基础环境配置gnina官方推荐在Linux环境下运行Ubuntu 20.04或者22.04是最省心的选择。我试过在CentOS 7上编译glibc版本太老折腾了很久还是放弃了。如果你用的是Windows建议直接上WSL2Ubuntu 22.04的镜像体验和原生Linux几乎没区别。macOS的话M1芯片的机器编译libmolgrid会有架构兼容问题Intel芯片的Mac倒是可以但性能不如Linux服务器。基础环境方面先确保你的系统有足够的磁盘空间。gnina编译过程中会产生不少中间文件加上CNN模型文件建议预留至少20GB。内存方面8GB是底线16GB以上会比较舒服因为CNN打分需要加载模型到内存里。GPU不是必须的但如果你有NVIDIA显卡编译CUDA版本会快很多。我实测过用CPU跑一个中等大小的蛋白-配体对接大概需要几分钟到十几分钟GPU的话可以缩短到几十秒。先更新系统包管理器然后安装编译工具链sudo apt update sudo apt install -y build-essential cmake git wget这几行命令看起来简单但有一个坑Ubuntu 22.04默认的gcc版本是11gnina的某些旧版本代码对gcc 11的支持不太好会报一些警告甚至错误。如果你遇到编译错误可以尝试安装gcc-9和g-9然后用update-alternatives切换过去。不过gnina的最新版本已经修复了这个问题所以优先用最新代码。2.2 libmolgrid的编译与安装libmolgrid是gnina的核心依赖之一它的作用是把分子结构转换成三维网格供CNN模型使用。这个库的安装是整个过程里最容易出问题的一步。官方文档给了编译步骤但有些细节没讲清楚。首先libmolgrid依赖Boost、OpenBabel和zlib。先装这些sudo apt install -y libboost-all-dev libopenbabel-dev zlib1g-dev注意Ubuntu仓库里的OpenBabel版本可能比较老gnina需要的是OpenBabel 3.0以上。如果你发现编译libmolgrid时报OpenBabel相关的错误可能需要从源码编译OpenBabel。源码编译OpenBabel的步骤稍微多一点但也不复杂wget https://github.com/openbabel/openbabel/archive/refs/tags/openbabel-3.1.1.tar.gz tar -xzf openbabel-3.1.1.tar.gz cd openbabel-openbabel-3.1.1 mkdir build cd build cmake .. -DCMAKE_INSTALL_PREFIX/usr/local make -j$(nproc) sudo make install编译完OpenBabel后需要设置环境变量让系统能找到它export LD_LIBRARY_PATH/usr/local/lib:$LD_LIBRARY_PATH export BABEL_LIBDIR/usr/local/lib/openbabel/3.1.0这两个环境变量很关键不设置的话后面编译libmolgrid时会报找不到OpenBabel库。建议把这两行加到~/.bashrc里省得每次都要手动设置。接下来编译libmolgridgit clone https://github.com/gnina/libmolgrid.git cd libmolgrid mkdir build cd build cmake .. -DCMAKE_INSTALL_PREFIX/usr/local -DOPENBABEL3_INCLUDE_DIR/usr/local/include/openbabel3 make -j$(nproc) sudo make install这里的-DOPENBABEL3_INCLUDE_DIR参数需要根据你的OpenBabel安装路径调整。如果你是用apt装的OpenBabel路径可能是/usr/include/openbabel3。cmake的时候注意看输出如果提示找不到OpenBabel就手动指定路径。注意libmolgrid编译过程中会下载一些测试数据如果你的网络环境不稳定可能会卡住。可以加上-DENABLE_TESTINGOFF来跳过测试数据的下载。2.3 CUDA与cuDNN的配置要点如果你打算用GPU加速CUDA和cuDNN的版本匹配是个大问题。gnina目前支持CUDA 10.2到11.8cuDNN需要对应版本。我建议用CUDA 11.8加cuDNN 8.6的组合这个组合在Ubuntu 22.04上最稳定。安装CUDA最省事的方法是用NVIDIA官方的runfile安装包注意安装时不要勾选驱动因为驱动最好单独用apt装。cuDNN的话下载对应的deb包用dpkg安装就行。安装完记得把CUDA的路径加到环境变量export PATH/usr/local/cuda-11.8/bin:$PATH export LD_LIBRARY_PATH/usr/local/cuda-11.8/lib64:$LD_LIBRARY_PATH验证CUDA是否装好可以用nvcc --version和nvidia-smi。如果nvidia-smi显示的CUDA版本和nvcc的不一致不用太担心只要驱动版本足够新就行。有一个坑我踩过gnina编译时默认会找系统里的CUDA但如果你装了多个版本的CUDAcmake可能会找错。可以在cmake时显式指定CUDA路径cmake .. -DCUDA_TOOLKIT_ROOT_DIR/usr/local/cuda-11.82.4 gnina主程序的编译与验证依赖都装好后编译gnina本身反而简单了git clone https://github.com/gnina/gnina.git cd gnina mkdir build cd build cmake .. -DCMAKE_INSTALL_PREFIX/usr/local -DCUDA_TOOLKIT_ROOT_DIR/usr/local/cuda-11.8 make -j$(nproc) sudo make install编译完成后运行gnina --help如果能看到帮助信息说明安装成功。如果报错说找不到libmolgrid.so检查一下/usr/local/lib是否在LD_LIBRARY_PATH里。可以用ldconfig来刷新动态库缓存sudo ldconfig还有一个常见问题gnina运行时需要CNN模型文件这些文件在编译时不会自动下载。你需要手动从gnina的GitHub release页面下载默认模型放到指定目录。默认模型文件是default2018.model和default2017.model放到gnina的安装目录下的models文件夹里或者用--cnn_model参数指定路径。3. 第一个对接案例从蛋白准备到结果分析3.1 蛋白与配体的预处理分子对接的第一步永远是准备受体和配体。gnina支持PDBQT格式的输入这和AutoDock Vina是一样的。如果你手头有PDB文件需要先转换成PDBQT。这里用OpenBabel或者AutoDockTools都可以。我习惯用OpenBabel的命令行工具速度快脚本化方便obabel receptor.pdb -O receptor.pdbqt -xr obabel ligand.sdf -O ligand.pdbqt -h-xr参数表示输出受体PDBQT-h表示加氢。注意OpenBabel加氢的规则和AutoDockTools不完全一样有时候会导致对接结果有差异。如果你对结果要求很严格建议用AutoDockTools的prepare_receptor4.py和prepare_ligand4.py脚本。蛋白预处理还有一个关键步骤去除水分子和无关的杂原子。gnina本身不会自动去水所以需要你手动处理。可以用grep命令过滤掉HOH开头的行grep -v HOH receptor.pdb receptor_nowater.pdb但要注意有些水分子是参与结合的盲目去掉可能会影响对接结果。我的经验是如果水分子在结合口袋内部并且和配体有氢键相互作用最好保留。可以用PyMOL或者Chimera可视化检查一下。配体预处理方面最重要的是生成正确的质子化状态和互变异构体。gnina的CNN打分对配体的质子化状态比较敏感如果pH不对打分会有偏差。可以用OpenBabel的--partialcharge参数来分配电荷或者用专门的工具如Epik、OpenEye的Omega来做质子化状态枚举。开源方案的话可以用Dimorphite-DL来生成不同pH下的质子化状态。3.2 定义对接盒子与参数设置对接盒子binding site的定义直接决定了对结果的好坏。gnina支持两种方式一种是直接指定盒子的中心坐标和尺寸另一种是用--autobox_ligand参数自动根据参考配体确定盒子。如果你有共晶结构里的参考配体用autobox最方便gnina -r receptor.pdbqt -l ligand.pdbqt --autobox_ligand reference_ligand.pdbqt --autobox_add 4--autobox_add 4表示在参考配体周围额外扩展4埃这个值可以根据口袋大小调整。如果口袋比较深可以加到6或者8。如果没有参考配体就需要手动指定盒子中心。可以用PyMOL或者Chimera找到口袋的中心坐标然后用--center_x、--center_y、--center_z指定用--size_x、--size_y、--size_z指定盒子尺寸。盒子尺寸一般设成20到25埃比较合适太小了配体可能跑不出去太大了计算量会增加。gnina的打分模式有几个选项最常用的是--cnn_scoring和--cnn_refinement。前者表示用CNN打分后者表示在Vina对接的基础上用CNN做精修。我一般用--cnn_scoring因为这样能充分发挥CNN的优势。如果你想要更快的速度可以用--cnn_refinement它先用Vina快速采样再用CNN对top构象重新打分。还有一个参数--exhaustiveness控制采样深度。默认是8对于虚拟筛选来说8已经够了。如果是精确对接可以调到16或者32但计算时间会成倍增加。我实测过exhaustiveness从8调到16对接时间大概增加一倍但RMSD小于2埃的构象比例提升不明显。所以除非你特别在意精度否则8就够用了。3.3 运行对接与结果解读参数设好后就可以运行对接了gnina -r receptor.pdbqt -l ligand.pdbqt --autobox_ligand reference_ligand.pdbqt --autobox_add 4 --cnn_scoring --exhaustiveness 8 --num_modes 9 --out output.sdf--num_modes 9表示输出9个构象这是Vina的默认值。gnina会输出一个SDF文件里面包含了每个构象的坐标和打分信息。SDF文件可以用PyMOL或者Chimera打开也可以用文本编辑器直接看。gnina的打分输出里有几个关键字段需要关注字段名含义解读要点CNNscoreCNN模型给出的打分范围0-1越高越好大于0.5通常认为有结合可能CNNaffinityCNN预测的结合亲和力单位是pIC50或pKd数值越大亲和力越强minimizedAffinityVina力场优化后的亲和力单位kcal/mol负值越大结合越强CNN_VS综合打分CNNscore和Vina打分的加权组合用于排序我一般用CNNscore来初筛用CNNaffinity来排序。CNNscore大于0.5的构象我会进一步用CNNaffinity来比较。如果两个构象的CNNscore差不多但CNNaffinity差很多那CNNaffinity高的那个更可信。结果分析时不要只看打分最高的那个构象。有时候打分第二、第三的构象在结合模式上更合理。我习惯把top 3的构象都导出来用PyMOL叠合一下看看它们的结合模式是否一致。如果top 3的构象结合模式差异很大说明这个配体的结合模式不确定需要谨慎对待。3.4 结果可视化与相互作用分析对接结果的可视化是验证对接质量的重要环节。我一般用PyMOL来做因为它脚本化方便可以批量处理。gnina输出的SDF文件可以直接用PyMOL打开pymol receptor.pdb output.sdf在PyMOL里可以用以下命令来显示相互作用# 显示蛋白口袋 show surface, receptor # 显示配体 show sticks, ligand # 显示氢键 distance hbonds, receptor, ligand, 3.5, mode2氢键的距离阈值一般设3.5埃这是供体和受体原子之间的距离。除了氢键还要看疏水相互作用、π-π堆积、盐桥等。这些相互作用可以用PLIP或者LigPlot来分析它们能自动生成相互作用图谱。有一个经验gnina的CNN打分对疏水相互作用的权重比较高所以如果配体主要是靠疏水作用结合的CNN打分通常会比较高。但如果配体主要靠氢键结合CNN打分可能不如Vina打分可靠。这时候可以结合两种打分来综合判断。4. 常见报错与性能调优的实战经验4.1 编译与运行时的典型报错gnina的安装和运行过程中有几个报错特别常见。我整理了一个速查表报错信息原因解决方法undefined reference toBabel::...OpenBabel链接错误检查LD_LIBRARY_PATH是否包含OpenBabel库路径CUDA error: no kernel image is availableCUDA架构不匹配在cmake时指定-DCUDA_ARCH_NAMEAlllibmolgrid.so: cannot open shared object file动态库路径未设置运行sudo ldconfig或手动export LD_LIBRARY_PATHCNN model file not found模型文件缺失从GitHub release下载模型文件用--cnn_model指定路径Segmentation fault during docking内存不足或输入文件格式错误检查PDBQT文件是否完整增加内存其中CUDA架构不匹配这个问题最隐蔽。gnina默认编译时会根据当前GPU的架构来生成代码但如果你把编译好的程序拿到另一台不同架构的GPU上跑就会报这个错。解决方法是在cmake时加上-DCUDA_ARCH_NAMEAll这样会为所有常见架构生成代码但编译时间会变长。还有一个坑gnina在运行时会尝试加载CNN模型如果模型文件路径不对它会静默失败然后用默认的Vina打分。这时候你看到的打分结果其实是Vina的不是CNN的。所以运行后一定要检查输出里有没有CNNscore字段如果没有说明CNN模型没加载成功。4.2 提升对接速度的实用技巧gnina的CNN打分比Vina慢这是事实。但通过一些技巧可以把速度提上来。第一个技巧是用--cnn_refinement模式。这个模式先用Vina快速采样只对top构象用CNN打分。实测下来速度比纯CNN打分快3到5倍而精度损失很小。对于大规模虚拟筛选这个模式是首选。第二个技巧是减少--num_modes。默认输出9个构象但如果你只关心最好的那个可以设成1或者3。这样能减少CNN打分的次数速度会快不少。第三个技巧是用GPU。gnina的CNN打分在GPU上比CPU快10倍以上。如果你有NVIDIA显卡一定要编译CUDA版本。我实测过一个包含1000个配体的库用CPU跑需要几个小时用GPU只要十几分钟。第四个技巧是并行化。gnina本身不支持多线程并行但你可以用GNU parallel或者xargs来同时跑多个配体。比如ls ligands/*.pdbqt | parallel -j 8 gnina -r receptor.pdbqt -l {} --autobox_ligand ref.pdbqt --cnn_scoring --out results/{/.}.sdf-j 8表示同时跑8个任务根据你的CPU核心数和GPU显存来调整。注意如果用的是GPU同时跑多个任务可能会导致显存不足需要根据显存大小来限制并行数。4.3 CNN打分与Vina打分的取舍策略gnina同时提供了CNN打分和Vina打分那到底该信哪个我的经验是分情况。对于虚拟筛选我倾向于用CNNscore来初筛因为CNNscore在区分活性化合物和非活性化合物方面通常比Vina打分好。但CNNscore的范围是0到1不同靶点的阈值不一样需要根据实际情况调整。我一般先用CNNscore大于0.5来筛一遍然后再用CNNaffinity来排序。对于结合模式预测我倾向于结合两种打分来看。如果CNN和Vina的打分都指向同一个构象那这个构象的可信度就很高。如果两者不一致就需要人工检查结合模式看看哪个更合理。还有一个细节gnina的CNN模型是在特定数据集上训练的对于某些靶点类别比如激酶、GPCR表现比较好。但对于金属酶、蛋白-蛋白相互作用界面表现可能不如传统打分函数。所以如果你的靶点比较特殊建议先用已知的共晶结构做一下验证看看gnina的打分是否合理。4.4 结果可重复性的保障措施分子对接的结果可重复性一直是个问题。gnina虽然比Vina好一些但也不是完全确定性的。为了保证结果可重复我一般会做几件事。第一固定随机种子。gnina的--seed参数可以指定随机种子默认是0。如果你想要完全可重复的结果可以设一个固定的种子比如--seed 42。但要注意即使种子固定不同硬件平台上的浮点运算差异也可能导致微小差别。第二记录完整的命令行参数。gnina的输出文件里会包含运行时的参数但有时候不全。我习惯把完整的命令行写到一个日志文件里方便以后复查。第三用多个随机种子跑几次看看结果是否一致。如果top构象的结合模式在多次运行中都很稳定那这个结果就比较可靠。如果每次跑出来的结合模式都不一样说明这个配体的结合模式不确定需要谨慎对待。第四保存中间文件。gnina在运行时会生成一些临时文件比如网格文件、模型文件等。这些文件对于复现结果很重要建议保留。5. 从单案例到批量筛选的扩展思路5.1 批量对接的脚本化实现跑通单个案例后下一步自然是批量对接。gnina本身没有提供批量处理的接口但用shell脚本或者Python脚本很容易实现。我一般用Python的subprocess模块来调用gnina这样可以灵活控制参数和并行度。下面是一个简单的批量对接脚本import os import subprocess from multiprocessing import Pool def dock(ligand): cmd [ gnina, -r, receptor.pdbqt, -l, ligand, --autobox_ligand, reference.pdbqt, --autobox_add, 4, --cnn_scoring, --exhaustiveness, 8, --num_modes, 3, --out, fresults/{os.path.basename(ligand)}.sdf ] subprocess.run(cmd, checkTrue) ligands [fligands/{f} for f in os.listdir(ligands) if f.endswith(.pdbqt)] with Pool(8) as p: p.map(dock, ligands)这个脚本会同时跑8个对接任务。注意如果用的是GPU需要根据显存大小调整Pool的大小。显存不够的话gnina会报CUDA out of memory错误。批量对接还有一个问题结果文件的组织。我习惯把每个配体的结果单独存一个SDF文件然后用一个汇总脚本提取打分信息生成一个CSV表格。这样方便后续用pandas做分析和排序。5.2 结果后处理与打分归一化批量对接完成后需要对结果进行后处理。gnina输出的SDF文件里包含了每个构象的打分但不同配体的打分范围可能不一样直接比较会有偏差。我一般会做几步归一化处理。第一步提取每个配体的top构象打分。可以用OpenBabel的obenergy或者自己写脚本解析SDF文件。我一般用RDKit来解析SDF提取CNNscore和CNNaffinity字段。第二步对CNNscore做z-score归一化。因为CNNscore的范围是0到1但不同靶点的分布不一样。z-score归一化后可以跨靶点比较。第三步对CNNaffinity做线性变换。CNNaffinity的单位是pIC50理论上可以直接比较但不同靶点的pIC50范围不一样。我一般会把CNNaffinity映射到0到1之间方便排序。第四步综合打分。我一般用CNNscore和CNNaffinity的加权平均作为最终排序依据权重可以根据实际情况调整。如果更看重结合可能性CNNscore的权重高一些如果更看重亲和力CNNaffinity的权重高一些。5.3 与商业软件的对比测试经验我用gnina和薛定谔的Glide做过对比测试用的是同一个靶点和同一批配体。结果发现gnina在早期富集上表现更好但Glide在结合模式预测上更准。这可能是因为Glide的打分函数经过了更精细的调参而gnina的CNN模型更擅长区分活性化合物和非活性化合物。具体来说在EF1%这个指标上gnina比Glide高了大概20%。但在RMSD小于2埃的构象比例上Glide比gnina高了大概10%。所以如果你的目标是虚拟筛选gnina是很好的选择如果你的目标是精确预测结合模式Glide可能更合适。当然这个对比只是基于一个靶点的测试不能代表所有情况。不同靶点、不同配体库结果可能不一样。我建议你在自己的靶点上做一下对比测试看看gnina是否适合你的需求。还有一个经验gnina的CNN模型对配体的分子量比较敏感。对于分子量小于300的配体CNN打分可能偏高对于分子量大于600的配体CNN打分可能偏低。所以如果你的配体库分子量分布很广建议按分子量分组来比较打分。5.4 实际项目中的部署建议如果你打算在课题组或者公司内部部署gnina有几个建议。第一用Docker或者Singularity来封装环境。gnina的依赖比较多手动安装容易出错。用容器可以保证环境一致性也方便迁移。gnina官方提供了Dockerfile可以直接用。第二建立标准化的输入输出流程。我一般会规定好PDBQT文件的命名规则、对接参数的默认值、结果文件的存放路径。这样不同人跑出来的结果可以互相比较。第三定期更新CNN模型。gnina的CNN模型会不定期更新新模型通常在更大的数据集上训练性能更好。建议每隔几个月检查一下GitHub上的release看看有没有新模型。第四做好日志记录。每次对接运行都要记录完整的命令行参数、运行时间、硬件配置、软件版本。这些信息对于复现结果和排查问题很重要。第五如果要做大规模虚拟筛选建议先用小规模测试集验证参数。我一般会选100个已知活性和100个已知非活性化合物跑一遍gnina看看富集因子是否合理。如果富集因子太低说明参数需要调整。6. 个人实操体会与后续扩展方向gnina的CNN打分版是我近几年用过的开源对接工具里最有意思的一个。它的出现打破了商业软件在高质量对接领域的垄断让预算有限的课题组也能做高质量的虚拟筛选。但gnina也不是万能的它的CNN模型有适用范围对于某些靶点类别可能不如传统打分函数。我在实际使用中发现gnina的CNN打分对结合口袋的疏水性比较敏感。如果口袋主要是疏水残基CNN打分通常比较准如果口袋有很多极性残基CNN打分可能偏高。这时候可以结合Vina打分来综合判断。还有一个体会gnina的--cnn_refinement模式在大多数情况下已经足够好了没必要非用--cnn_scoring。--cnn_refinement的速度快很多精度损失很小特别适合大规模筛选。只有在需要最高精度的时候才用--cnn_scoring。后续扩展方向的话我觉得有几个值得尝试。一是把gnina的CNN打分和分子动力学模拟结合起来用MD来验证对接结果的稳定性。二是把gnina集成到自动化流程里比如用Snakemake或者Nextflow来管理整个虚拟筛选流程。三是尝试训练自己的CNN模型用自己靶点的数据来微调gnina的默认模型可能会得到更好的结果。最后分享一个小技巧gnina的输出SDF文件里每个构象的CNNscore和CNNaffinity都在注释行里。如果你用PyMOL打开SDF可以用iterate命令来提取这些打分iterate all, print(name, resn, CNNscore, CNNaffinity)这样可以快速比较不同构象的打分不用一个个手动看。