1. 先搞清楚DOCK6是什么再谈怎么用分子对接这个词搞计算化学和药物设计的人肯定不陌生。简单说就是把一个小分子配体放到一个大分子受体的结合位点里通过几何匹配和能量打分预测它们最可能的结合模式和亲和力。UCSF DOCK6就是干这个事的经典开源程序之一来自加州大学旧金山分校的Irwin和Shoichet团队从九十年代一路迭代到现在目前最新的6.x版本依然活跃在各类虚拟筛选和药物先导发现的项目里。很多人一开始会纠结现在AutoDock Vina、Glide、GOLD都很好用为什么还要学DOCK6我的答案是DOCK6的定位和它们不太一样。DOCK6的底层逻辑是“锚定搜索”anchor-and-grow它先把配体拆成刚性锚用几何匹配的方式把锚放到口袋中再逐步长出柔性部分。这套逻辑在应对构象搜索和结合模式多样性时非常稳尤其在处理多构象、大分子库的虚拟筛选时预计算好的网格配合上高效打分速度快、可解释性强而且完全免费、开源、无需License对于学生党和小课题组特别友好。这篇教程面向的读者很明确第一次接触DOCK6的人或者之前只用过Vina、想补充一套更严谨的对接流程的人。我会按照实际项目的完整路径来讲——从安装、受体准备、球集生成、网格计算、配体准备到dock.in编写、运行对接、结果分析最后把我在实操中踩过的坑全部列出来。有的细节官方文档语焉不详网上教程又互相矛盾我尽量把“为什么这么做”也一并讲清楚而不是只丢给你一串命令。2. 环境准备与安装DOCK6编译其实没那么难2.1 依赖安装一次到位DOCK6的安装不算复杂但因为它需要自己编译很多新手会卡在依赖上。在Linux环境Ubuntu/CentOS都行下建议先把这几样装齐gcc/g、make、flex、bison、zlib开发包。如果你还想用GPU加速版本那还要准备CUDA Toolkit和对应的显卡驱动。Ubuntu/Debian系统一行命令就能搞定大部分依赖sudo apt update sudo apt install build-essential flex bison zlib1g-dev如果你的系统是CentOS/RHEL对应的是sudo yum install gcc gcc-c make flex bison zlib-devel这里多说一句flex和bison这两个缺一不可它们是用来解析DOCK6输入文件的如果没装编译到一半会报一堆parser错误别问我怎么知道的我第一次编译时就漏了bison卡了整整一个下午。2.2 编译与自检去UCSF DOCK官网下载对应版本的源码压缩包解压后进入目录按顺序执行tar -xzvf dock6.tar.gz cd dock6 ./configure gmake gmake test./configure会检测你的编译环境看到“Configuration successful”之类的提示再继续。gmake就是正式编译这个过程根据机器配置可能要几分钟到十几分钟属于正常现象。编译结束后gmake test会跑官方自带的小测试集这一步强烈建议不要跳过它能帮你确认打分函数、网格计算这些模块是否都正常。如果测试全部通过你会看到类似“test completed successfully”的输出。此时可以在dock6目录下找到bin/dock6、bin/grid、bin/sphere_generator、bin/dms、bin/showbox等一系列可执行文件建议把它们所在的bin目录加进PATH后面操作会省很多事。2.3 自测用例花十分钟验证环境编译完成后最快验证安装是否可用是直接在解压目录下找一个官方测试例跑一遍。比如DOCK6自带test/目录里面有针对不同功能的测试用例随便挑一个cd进去执行../../bin/dock6 -i dock.in -o dock.out看dock.out最后几行有没有正常的打分输出或者直接查看输出的分子文件是否生成。这一步如果能顺利跑通说明你的软件环境基本没问题后面出问题就可以甩锅给结构文件而不是编译环境了。3. 受体准备所有坑的大本营我见过太多人花大力气把DOCK6装好结果对接出来的东西完全不能用问题就出在受体准备上。DOCK6对受体结构的要求比Vina严因为它要拿受体来做网格打分原子类型、电荷、加氢状态都会直接影响静电项和范德华项。3.1 获取并清洗PDB结构从RCSB PDB下载晶体结构后第一步不是急着加氢而是先把体系清理干净。晶体结构里通常含有水分子、配体、金属离子、糖链、去污剂等哪些该保留哪些该删取决于你的研究目标。常规做法是只保留受体蛋白本身grep ^ATOM receptor.pdb receptor_protein.pdb如果你结合位点里有一个结构水分子对氢键至关重要也可以手动把它加回来。金属离子如果是功能必需的比如锌依赖的金属蛋白酶一般建议保留并在后续电荷分配时做特殊处理。但如果你是新手最简单的策略是先全部删掉水只保留蛋白链等基本流程跑通了再考虑细节。晶体结构中常缺失原子或残基如果缺失部分正好在结合位点附近建议先去补环可以用ChimeraX的Model Loop或者Modeller否则对接结果基本不能用。这个前期检查值得做省得后面白跑。3.2 加氢和电荷一步都不能省PDB结构几乎都不含氢原子而DOCK6做分子力学打分时不能没有氢所以加氢是硬性要求。推荐用DOCK6自带的reduce工具或者用ChimeraX的Dock Prep、PDB2PQR都行。我的习惯是用ChimeraX做加氢和局部矫正因为可视化方便能顺便检查His的翻转、Asn/Gln侧链是否需要旋转。加完氢之后最关键的一步给受体分配电荷。这一步很多人会漏后果就是静电打分整个就是错的。DOCK6官方教程比较推荐的做法是先把含氢结构转成mol2文件再借助AMBER的antechamber工具分配AM1-BCC电荷或Gasteiger电荷。在ChimeraX里可以这样操作打开加氢后的结构Tools → Structure Editing → Dock Prep勾选“Add charges”和“Write mol2”保存时选择AMBER ff14SB力场这样输出的mol2文件带上了正式的原子类型和电荷。如果你偏好纯命令行也可以用antechamberantechamber -i receptor_h.pdb -fi pdb -o receptor.mol2 -fo mol2 -c bcc -nc 0 -s 2注意-nc 0表示受体净电荷为0如果你知道蛋白实际带净电荷要对应修改。另外DOCK6在做网格计算时需要的是同时包含PDB和mol2两套文件PDB用于生成表面和球集mol2用于网格电荷。3.3 加氢里暗藏的坑加氢这件事看着简单其实有很多细节。比如His残基有三种质子化状态HID、HIE、HIPpH环境不同结合位点里His的质子化状态差异会直接影响氢键网络和打分。晶体结构本身无法区分这几种状态需要你根据周围环境和目标pH去判断。还有一个我吃过亏的地方ChimeraX加氢默认会把所有可电离残基按生理pH处理但有些结合口袋pH环境很特殊比如某些溶酶体酶作用环境偏酸如果你不加思考直接用默认状态静电打分会偏差很大。所以做之前最好查阅文献确认结合位点的关键残基质子化状态。4. 生成分子表面与球集DOCK6的几何搜索核心这一章节是DOCK6区别于其他对接软件的核心也是新手最容易一头雾水的地方。DOCK6不像Vina那样用盒子定义搜索空间然后全局搜索它用“球集几何匹配”的方式把配体放到口袋里所以球集的质量直接决定对接成败。4.1 dms先生成分子表面球集的生成需要两个前置信息受体分子的溶剂可及表面以及这个表面上每个点的法向量。DOCK6调用dms程序来做这件事。命令大概是这样的$DOCK6/bin/dms receptor_protein.pdb -n -w 1.4 -v -o receptor.dms解释一下参数-n表示忽略HETATM也就是忽略配体、水等非蛋白原子如果你把关键水或离子保留下来了这个参数要不要加得慎重-w 1.4是探针半径模拟水分子的大小-v是输出详细信息-o指定输出文件名。生成的receptor.dms是一个点云文件记录了表面上每个点的坐标和法向。如果这一步报错最常见的原因是PDB文件里原子类型不规范或者结构里有原子坐标异常。建议先检查PDB是否只含有标准氨基酸非标准残基或错误元素符号都会让dms崩溃。4.2 sphere_generator从点云到球集有了表面点云下一步就是用sphere_generator生成“负镜像球”。这个程序会沿着表面点法向量朝口袋内部放置球体球的大小和位置描述了结合腔的几何形状。$DOCK6/bin/sphere_generator -i receptor.dms -o receptor.sph -a 640 -p 0.0 -r 1.4 -l 5参数含义大概是这样-a 640是单次簇球数量上限-r 1.4是探针半径-l 5是层数限制-p是偏移距离。不同版本的默认值略有差异但上述这套参数是很多教程通行的用法可以直接用。生成的receptor.sph就是包含所有可能球体的文件但里面球太多了必须筛选。这里就体现出参考配体的价值如果你对接的对象已有共晶配体可以用sphere_selector把结合位点附近的球挑出来$DOCK6/bin/sphere_selector receptor.sph ligand_original.mol2 10.0这条命令会保留距离参考配体10 Å以内的球体输出文件默认叫selected_spheres.sph。如果你没有参考配体只有已知的活性口袋残基也可以手动计算口袋中心再用类似的距离阈值做筛选。4.3 球簇选择的实操经验筛选球集时最直接的经验是球集不要贪多。官方教程和很多案例里选取15到30个球就够了覆盖口袋的主要子区域即可。球太多会导致后续几何匹配时搜索空间暴涨速度变慢而且容易产生假阳性结合模式。用ChimeraX或PyMOL把selected_spheres.sph和受体叠在一起看是查球集质量最直观的方式。合格的球集应该像一层“假想配体骨架”贴合在结合腔里疏水口袋里有球覆盖极性区域也有球但相对稀疏。如果你发现某些球穿出了蛋白表面或者离口袋十万八千里说明dms表面计算有误或者选择阈值设得太大。我个人的经验是宁可多花10分钟在PyMOL里审视球集也不要直接带着一个有问题的球集去跑对接因为后面所有的结果都会被带偏而且很难排查。5. 网格计算为打分函数铺路5.1 盒子怎么划DOCK6的打分依赖预计算的网格范德华网格和静电网格。网格覆盖的区域必须包含整个结合口袋同时不能太大否则计算量陡增且打分项失去局部性。生成盒子最标准的做法是用showbox程序。你需要手动创建一个box.in文件内容类似YES receptor_ligand.pdb box.pdb 20 20 20第一行YES表示用配体计算盒子中心第二行是包含蛋白和参考配体的结构文件第三行是输出盒子文件的名称最后三行是盒子在x、y、z三个方向上的半边长单位Å。然后运行$DOCK6/bin/showbox box.in生成的box.pdb就是网格盒子可以在ChimeraX里和受体、配体叠放检查。常规经验是盒子边长要能覆盖配体四周至少3到5 Å的余量对接小分子的话半边20 Å通常够用但如果你筛的是多肽或大分子片段要相应扩大。5.2 grid.in参数拆解网格计算用grid程序输入文件叫grid.in。一个常用的模板如下compute_grids yes grid_spacing 0.4 output_molecule no contact_score yes contact_cutoff_distance 4.5 energy_score yes energy_cutoff_distance 999 atom_model a attractive_exponent 6 repulsive_exponent 12 distance_dielectric yes dielectric_factor 4 bump_filter yes bump_grid_prefix bump receptor_out_file receptor_box.pdb box_file box.pdb vdw_grid_file vdw es_grid_file es其中有几个参数值得重点说说。grid_spacing是网格间距单位Å。0.4 Å是推荐值兼顾速度和精度。如果你的体系很大、只想粗筛可以放宽到0.5如果做精细的相互作用分析可以试试0.3但内存和时间成本会明显上升。attractive_exponent和repulsive_exponent控制范德华势函数的形状6-12是经典Lennard-Jones形式这也是AMBER力场的标配。distance_dielectric yes加上dielectric_factor 4表示静电用距离相关介电常数这是很多蛋白-配体打分里降低长程静电影响的常见做法。bump_filter yes会额外生成一个用于碰撞检测的网格这个网格在对接时用来快速排除那些和受体有严重空间冲突的配体构象。vdw_grid_file和es_grid_file是输出网格文件的前缀之后对接程序会自动读取。运行命令$DOCK6/bin/grid -i grid.in -o grid.out成功后会生成vdw.bmp、es.bmp等网格文件以及一个带盒子的受体参考PDB。5.3 网格质量的检查方法网格算完至少要做三件事第一确认grid.out里没有报错或警告比如“atom type unknown”这类第二检查receptor_box.pdb的盒子是否完整覆盖了球集区域第三有条件的话把参考配体放到盒子里看它的能量是否表现为负值如果参考配体在结合位点里能量还是正的多半是电荷或原子类型分配出了问题。我遇到过最诡异的情况是网格文件生成正常但跑到对接时所有配体打分都异常高查了半天发现是bump_grid_prefix里的bump网格没算出来导致后续碰撞过滤逻辑出错。所以网格计算这块仔细看输出日志永远不会错。6. 配体准备从二维结构到三维构象配体的准备和受体同样重要但经常被新手忽略。很多人直接从PubChem下载一个SDF就往DOCK6里塞结果各种报错或者打分离谱。配体的准备流程包括三维结构生成、加氢和电荷、构象搜索三步。6.1 生成三维结构如果你手里只有SMILES或者二维结构先用Open Babel或RDKit生成三维坐标obabel ligand.smi -O ligand_3d.mol2 --gen3dRDKit用户在Python里用ETKDG方法也行from rdkit import Chem from rdkit.Chem import AllChem mol Chem.MolFromSmiles(你的SMILES) mol Chem.AddHs(mol) AllChem.EmbedMolecule(mol, AllChem.ETKDG())注意检查生成的三维结构有没有不合理的环张力或者原子重叠必要时用MMFF或UFF做一轮能量最小化。6.2 加氢与电荷分配DOCK6要求配体mol2文件里有明确的SYBYL原子类型和电荷。最通用的做法是用AMBER的antechamberantechamber -i ligand_3d.mol2 -fi mol2 -o ligand.mol2 -fo mol2 -c bcc -nc 0 -s 2-c bcc用AM1-BCC电荷对药物类小分子比较合适-nc根据配体的质子化形式设置净电荷比如羧基去质子化后就是-1。在生理pH下该用哪个质子化状态强烈建议先用DataWarrior或ChemAxon做个预测不要想当然地认为中性分子就是不带电的。6.3 构象搜索DOCK6对接前最值得花时间的环节分子对接的核心难点之一是构象搜索。即使DOCK6有anchor-and-grow的柔性搜索能力配体的起始构象仍然会影响最终结果尤其是那些环系多、柔性大的分子。我建议在对接前用RDKit或OMEGA生成一组多构象集合然后把它们合并成一个多分子的mol2或SDF文件再喂给DOCK6。RDKit的做法很直接from rdkit import Chem from rdkit.Chem import AllChem params AllChem.ETKDGv3() params.randomSeed 42 mol Chem.AddHs(Chem.MolFromSmiles(你的SMILES)) AllChem.EmbedMultipleConfs(mol, numConfs100, paramsparams)也可以直接用Open Babel的confab模式obabel ligand.mol2 -O ligand_confs.sdf --conformer --nconf 100 --score rmsd --rmsd 0.5--rmsd 0.5表示构象之间RMSD小于0.5 Å的视为重复并剔除避免输出一堆几乎一模一样的冗余构象。这一步生成的多样构象会显著提高DOCK6找到合理结合模式的概率。7. 编写dock.in并正式对接7.1 一个能直接跑的柔性对接模板DOCK6的输入文件dock.in用的是“关键字 值”的写法没有复杂的标记语言但这个文件踩的坑也不少。下面是一个经过验证的柔性对接模板以网格打分为主含能量最小化ligand_atom_file ligand_confs.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd yes use_database_filter no orient_ligand yes automated_matching yes receptor_site_file selected_spheres.sph max_orientations 1000 critical_points 1 chemical_matching no use_ligand_spheres no bump_filter yes score_molecules yes contact_score_primary no contact_score_secondary no grid_score_primary yes grid_score_secondary no grid_score_rep_constant 1 grid_score_lin_constant 1 grid_score_vdw_scale 1 grid_score_es_scale 1 grid_score_att_exp 6 grid_score_rep_exp 12 minimize_ligand yes minimize_anchor yes minimize_flexible_grow yes use_advanced_minimizer yes use_initial_energy_only yes minimization_algorithm simplex simplex_max_iterations 1000 simplex_tors_premin_iterations 20 simplex_max_grow_iterations 500 simplex_initial_convergence 1.0 simplex_convergence 0.1 simplex_score_converge yes simplex_cycle_max 1 output_molecule yes output_mode all output_filename docked_ligands.mol2这个模板把conformer_search_type对应的柔性搜索、锚定增长、能量最小化都打开了。注意我没写conformer_search_type这一行是因为DOCK6默认就会根据orient_ligand和automated_matching的组合进入柔性搜索模式如果你读到的DOCK6版本提示找不到这个参数不用慌上述字段对齐即可。7.2 关键参数解释别乱调的清单receptor_site_file指向你筛选后的球集文件这个必须和受体对应否则对接会直接在错误的位置搜索。max_orientations是配体锚在结合位点里的最大取向数默认1000一般够用。你要是发现结果里构象总是只有那几种可以适当提高到2000到5000但也要明白取向越多计算越慢收益有时并不是线性的。critical_points表示球集中至少要匹配几个点才算有效取向。设为1最灵活取2或3会更严格但有可能漏掉好的结合模式。对新手我推荐先设1后续根据结果再调整。打分方面grid_score_primary yes表示最终排序以网格打分为主也就是范德华和静电的加和。contact_score只是辅助用于评估形状互补性。minimize_ligand、minimize_anchor、minimize_flexible_grow这三个都建议打开这会在粗打分之后对配体做一轮局部能量优化明显减少立体冲突和键长键角不合理的情况。simplex_*一系列参数控制最小化算法。DOCK6内置了简单下坡法虽然不是最精细的优化器但在对接场景下够用而且稳定。simplex_score_converge yes表示能量变化低于收敛阈值就提前停能省不少时间。7.3 运行与输出命令非常简单$DOCK6/bin/dock6 -i dock.in -o dock.out跑的过程中可以随时打开dock.out看进度它会逐步输出每个分子的取向匹配、打分、最小化信息。结束后在docked_ligands.mol2文件里能看到所有保留的对接构象每个分子带一个Grid_Score或类似名称的注释就是网格打分值分数越低表示结合越有利。如果你准备跑虚拟筛选可以把多个配体分子合并成一个多分子mol2或SDF文件dock.in里的ligand_atom_file指向这个集合文件即可。limit_max_ligands no表示不限制处理数量用于筛选上万分子的数据库。8. 结果分析RMSD、打分与构象聚类8.1 先判断对接是否成功拿到dock.out之后第一件事是看有没有Error和Warning。常见的问题是“no orientations found for molecule X”意思是某个配体在球集附近找不到合法的取向大部分时候是因为球集选得不好而不是配体本身的问题。再看docked_ligands.mol2里的打分分布。如果所有打分都是很大的正值先别急着做结论回查网格和电荷如果分数集中在负值到零附近说明整体环境是合理的可以继续后面的构象和后处理分析。8.2 RMSD有共晶结构的人必查如果你手里有参考共晶配体那么calculate_rmsd yes会在每个输出分子的注释里记录它对参考配体的RMSD值。挑RMSD最低的打分构象如果RMSD小于2 Å说明你的流程和方法基本到位如果所有低打分构象RMSD都大于4 Å那就要警惕了——可能是球集选偏、受体质子化状态错误也可能是配体初始构象和实际结合构象差异太大。这里有个小技巧用PyMOL把对接后的最优pose和共晶配体叠合肉眼看相互作用模式氢键、疏水接触。RMSD数字只能说明原子坐标差异而氢键网络是否合理、疏水残基是否被利用才是评价结合模式是否有说服力的关键。8.3 构象聚类和可视化DOCK6输出的构象多的时候可能上百个直接一个个看会崩溃。建议用构象聚类工具按RMSD聚类比如ChimeraX自带的相关分析或者用Python脚本处理。一般取能量最低、且在不同聚类里代表性能最好的那几个pose作为后续分子动力学模拟或结合自由能计算的起点。可视化的常规套路是在ChimeraX里打开受体结构把球集、对接pose、关键残基的侧链都显示出来用氢键距离标注好导出高清图。论文里的对接图够用的标准是口袋表面光滑、配体清爽、关键相互作用有距离标线。9. 常见问题与避坑速查表以下是我在实际使用DOCK6过程中遇到过的典型问题整理成速查表建议收藏。现象可能原因解决办法编译报错找不到flex/bison依赖未安装安装flex和bison后重新./configure和gmakedms程序崩溃或输出为空PDB含有非标准原子或坐标异常清洗PDB只保留标准氨基酸和必要辅因子生成球集数量异常多且乱dms表面法向计算错误检查dms参数中的探针半径用可视化确认表面质量球集不在结合位点参考配体坐标未和受体对齐盒子中心错误检查参考配体位置使用sphere_selector时保证配体和受体的坐标体系一致grid报atom type unknownmol2原子类型不规范用antechamber或ChimeraX重新生成mol2确保SYBYL原子类型完整对接时no orientations found球集太稀疏或配体太大增加球集数量检查球集是否覆盖整个配体适当增大max_orientations所有配体打分都是正数电荷分配错误或网格盒子未覆盖口袋检查配体/受体电荷重新生成网格扩大盒子范围输出mol2文件为空dock.in中output_molecule未开启或write_frequency设置过大将output_molecule设为yes输出模式设为all对接结果高分的RMSD反而大打分函数局限或球集引导了错误取向调整球集筛选增加构象搜索和最小化尝试不同初始构象虚拟筛选速度太慢网格间距过密或取向数过大网格间距改为0.5适当降低max_orientations如果你用的是较新版本的DOCK66.9及以上部分程序的参数名可能和我上面给的略有出入以官方release note和测试目录里的输入文件为准。遇到参数解析失败最简单的方法是把官方test目录里对应功能的dock.in或grid.in打开对照着改。最后再分享一个小技巧DOCK6是一个高度模块化的程序新手最容易把它当成一个“黑盒按钮”来用但它真正的价值恰恰在于每一步都是透明的——球集、网格、打分都可以单独检查和可视化。我每次跑新项目都会强制自己把球集和网格的可视化检查纳入流程虽然多花十分钟但换来的是后续少熬几个夜。这套流程你多走几遍之后会发现它比那些“一键对接”的工具更能帮你建立对分子识别的直觉。