拿到一个新的靶标蛋白,第一反应不是急着找抑制剂,而是先问自己一句:这个蛋白到底在哪儿结合小分子?口袋有多大?氨基酸残基长什么样?这些问题靠猜没用,分子对接就是用来回答这件事的。AutoDock Vina 是目前分子对接领域最常用的开源程序之一,安装简单、参数不复杂、单机也能跑,很多药物设计、虚拟筛选的入门项目都用它起步。
这篇文章我会以 PDB 数据库里的 7VU6 结构为案例,带你完整跑一遍 AutoDock Vina 的分子对接流程。从下载结构、预处理受体和配体,到设置搜索盒子、运行 Vina、解读打分结果,每一步都会拆开讲,顺便把我自己踩过的一些坑写出来。不管你是刚接触计算化学的学生,还是想给实验室课题加一个计算验证的科研人员,这套流程都可以直接照抄。
1. 先把AutoDock Vina这套工具讲透
1.1 分子对接到底在解决什么问题
分子对接,简单说就是预测一个小分子配体和一个大分子受体之间最可能的结合方式。蛋白好比一把锁,小分子好比钥匙,对接算法就是在受力平衡、空间互补、化学匹配这些约束下,搜索钥匙插入锁孔的最优姿态。这个“结合方式”不只是位置和方向,还包括小分子的柔性构象变化,也就是可能出现的不同姿势。
AutoDock Vina 采用的是基于经验的打分函数,它把范德华相互作用、氢键、静电、去溶剂化、熵效应等因素通过回归拟合转化为一个打分值,单位是 kcal/mol。这个分数越负,理论上结合越有利。Vina 做得好的地方在于搜索速度快、打分精度在同等级工具里性价比非常高,它不像分子动力学那样需要几十微秒的轨迹,也不像对接刚起步时的软件那样动不动就卡死。你只需要给出受体结构、配体结构和搜索空间,它就能在几分钟到几个小时内给你一批候选构象。
当然,对接不是终点。它解决的是“能怎么结合”的问题,至于“结合得多牢”“会不会在生理条件下真正发生”,还需要后续的分子动力学、自由能计算甚至湿实验去验证。但 Vina 作为第一层过滤器,完全够用。
1.2 为什么是AutoDock Vina
目前市面上的对接工具不少,比如 Glide、GOLD、MOE,这些商业软件精度高但价格高,而且许可证管理麻烦。AutoDock Vina 的优势有三个:
第一,免费开源。对个人学者和小课题组来说,不用写预算本,下载就能用。第二,跨平台支持好。Windows、Linux、macOS 都有编译好的版本,命令行工具对脚本化批量操作非常友好。第三,社区庞大。网上教程多,PDB 数据、打分函数、参数选择的各种经验,很容易搜到。
和它的前身 AutoDock 4 相比,Vina 在搜索算法和打分函数上都做了大幅简化。AutoDock 4 需要用户手动指定许多原子类型和 torsion tree,而 Vina 对大多数标准有机小分子只要把 PDBQT 文件准备好就行。它还有一个很实用的特点:会自动处理配体的可旋转键,不需要你手动定义柔性残基,除非你做的是柔性对接。这让入门门槛低了很多。
但我必须提醒一句:Vina 不是万能药。它对高度柔性的配体、金属蛋白、共价抑制剂这类体系,预测能力有限。尤其是含金属离子的靶标,Vina 的打分函数没有很好地处理配位键,结果可能偏差很大。如果遇到这种体系,要么换工具,要么用专门的金属蛋白对接方案。
1.3 一份能直接照抄的流程概览
做 AutoDock Vina 对接,完整流程可以拆成六个环节:
- 获取受体结构文件,一般从 PDB 数据库下载。
- 从复合物结构中提取配体,或者自己画配体并生成 3D 构象。
- 对受体、配体分别做预处理,加氢、去水、合并非极性氢、加电荷、转成 PDBQT 格式。
- 确定对接盒子:中心坐标和边长尺寸。
- 运行 Vina,设置搜索精度(exhaustiveness)和输出模式数。
- 分析结果:看打分值、聚类情况、氢键和疏水接触,再用 PyMOL 或 ChimeraX 可视化。
这套流程里,最容易被忽略也最影响结果的是第 4 步。盒子位置错了,后面所有对接都是废的。很多新手跑出来的分数很漂亮,但姿态明显不合常理,十有八九就是盒子没放对。
2. 环境搭建与文件准备
2.1 安装AutoDock Vina和配套工具
AutoDock Vina 的安装没有那么玄学。Linux 环境下最省事的方式是从官网下载编译好的二进制包,解压后直接运行。
wget https://github.com/ccsb-scripps/AutoDock-Vina/releases/download/v1.2.5/vina_1.2.5_linux_x86_64 chmod +x vina_1.2.5_linux_x86_64 sudo mv vina_1.2.5_linux_x86_64 /usr/local/bin/vina vina --versionWindows 用户可以从同样地方下载 exe 文件,放到一个不含空格的路径下,用 PowerShell 或者 CMD 调用。macOS 的 M 系列芯片需要注意选择对应的 arm64 版本,别下载成 x86_64,虽然 Rosetta 也能跑,但原生版本更快。
除了 Vina 本体,你还需要准备两套配套工具。
第一套是 MGLTools 里的 AutoDockTools,主要用来处理受体和配体的 PDBQT 文件。MGLTools 的安装包里包含可以独立调用的 Python 脚本,比如 prepare_receptor4.py 和 prepare_ligand4.py。新版 MGLTools 已经支持 Python 3,直接到官网下载对应平台的安装包即可。
第二套是 Open Babel,它的作用非常广,格式转换、加氢、计算 Gasteiger 电荷都能做。很多人觉得 PDBQT 文件不好写,Open Babel 一行命令就能从 mol2、SDF、SMILES 转出 PDBQT。在 Ubuntu 上安装是:
sudo apt install openbabel也可以从 conda 安装:
conda install -c conda-forge openbabel如果后续做批量虚拟筛选,建议把 Vina 和 Open Babel 都配到系统 PATH 里,这样写循环脚本会舒服很多。
2.2 下载7VU6结构并做初步检查
PDB 数据库的每个结构都有一个四位字符 ID,7VU6 就是其中一个以 7 开头的较新条目。打开 RCSB PDB 网站,搜索 7VU6,在 Structure 页面下载 PDB 格式文件,文件名一般叫 7vu6.pdb。
下载之后,不要直接拿来对接,先做三件事。
第一,用文本编辑器打开 PDB 文件,看这个结构是不是复合物。如果里面既有蛋白链,又有 HETATM 开头的小分子,说明包含了共晶配体,这种结构做对接再合适不过。你可以用共晶配体的坐标来确定结合口袋。如果只有蛋白没有配体,就需要通过文献或预测工具来找口袋,工作量大一些。
第二,看分辨率。晶体结构的分辨率是判断结构质量的重要指标。分辨率小于 2.0 的蛋白结构基本靠谱,2.0 到 2.5 要看电子密度图情况,大于 3.0 的结构对接时就得格外小心,侧链位置的误差可能直接导致口袋变形。7VU6 如果是近两年解析的结构,一般分辨率不会太差,但还是要自己确认。
第三,看有没有缺失残基。很多 PDB 结构因为柔性区域太乱,无法解析出完整的电子密度,文件里会有 REMARK 465 记录缺失残基。缺失的位置如果在口袋附近,说明该结构对结合位点建模不完整,可能不适合直接对接。遇到这种情况,最好用同源建模补全,或者换一个更高质量的结构。
下载完成后,建议用 PyMOL 或 ChimeraX 快速看一眼 7VU6 的整体形貌,确认蛋白链数量、配体位置、水分子数量。这个步骤虽然简单,却能在后续预处理时节省很多时间。
2.3 受体与配体文件格式转换
AutoDock Vina 要求的输入格式不是普通 PDB,而是 PDBQT。它是在 PDB 基础上扩展出来的格式,多了两列关键信息:原子类型和 Gasteiger 电荷。Vina 的打分函数依赖这些信息来计算相互作用,所以转换这一步不能省,也不能乱省。
PDBQT 里每个原子行的末尾会有类似A C这样的原子类型标记。碳原子会被识别为 A(芳香碳)、C(脂肪碳)等,氮、氧、硫、磷也各有对应类型。如果文件里出现 Vina 不认识的原子类型,程序会直接报错。这也是为什么预处理工具不能随意跳过。
转换受体文件最常规的命令是:
python prepare_receptor4.py -r 7vu6.pdb -A hydrogens -U nphs_lps_waters -o receptor.pdbqt这个命令里的-A hydrogens表示给受体加氢。-U nphs_lps_waters表示移除非极性氢、孤对电子和水分子。处理后输出的 receptor.pdbqt 就是干净可用的受体文件。
配体转换可以用 Open Babel:
obabel ligand.mol2 -O ligand.pdbqt --partialcharge gasteiger注意,Open Babel 转换 PDBQT 时对某些原子的类型判定可能不完美,比如有机磷、磺酸基这类特殊基团,容易转出 Vina 不认识的原子类型。如果你发现对接时总报Unknown atom type的错,可以改用 MGLTools 的 prepare_ligand4.py:
python prepare_ligand4.py -l ligand.mol2 -o ligand.pdbqt这个脚本会按照 AutoDock 的原子类型规则重新标记,处理更靠谱。
3. 四个关键参数,定下对接的边界
3.1 盒子坐标:中心在哪、尺寸多大
参数里最要命的就是盒子(search space)。Vina 只会在你指定的长方体范围内搜索配体的结合位置和姿态。如果盒子不包含真实结合位点,哪怕打分再低也没意义,因为那个姿势分子根本不可能在蛋白里存在。
确定盒子中心,首选方法是使用共晶配体的中心。7VU6 里如果带有配体,直接在 PDB 文件里找 HETATM 记录,计算所有配体原子的坐标平均值,这个质心就是最合理的盒子中心。可以用一段简单的 Python 脚本算:
import sys coords = [] for line in open(sys.argv[1]): if line.startswith(("ATOM", "HETATM")): coords.append((float(line[30:38]), float(line[38:46]), float(line[46:54]))) n = len(coords) cx = sum(c[0] for c in coords) / n cy = sum(c[1] for c in coords) / n cz = sum(c[2] for c in coords) / n print(f"{cx:.3f} {cy:.3f} {cz:.3f}")把配体坐标文件传给脚本,拿到中心坐标。接着设置盒子大小,一般用size_x、size_y、size_z三个边长。对于单配体结合口袋,边长 20 埃左右比较合适。如果结合位点结构未知,只能做盲对接,尺寸就需要扩大到覆盖整个蛋白,甚至 60 到 80 埃,但这样搜索空间大,时间会成倍增长,结果也不一定稳定。
我个人的经验是:一旦你有了共晶配体,就把盒子尺寸控制在刚好包住配体外加 5 到 8 埃的余量,这样可以避免 Vina 把配体推到口袋外面,也能减少不合理的远距离构象。
3.2 exhaustiveness:搜索力度与耗时的折中
exhaustiveness 是 Vina 最重要的搜索参数,它控制全局搜索的随机尝试次数。数值越高,搜索越充分,但耗时也越长。默认值是 8,适合快速测试。我的建议是正式对接至少用 16,条件允许用 32,如果要发文章或者做关键预测,直接上 64。
有人可能会问,是不是越高越好?理论上是的,但收益会递减。对一个 20 埃见方的盒子,exhaustiveness 从 8 提升到 32,往往能找到明显更优的构象;但从 32 提升到 64,差别就不大了,耗时却翻倍。所以在批量筛选时,我通常先用 8 跑一遍快速过滤,把明显不合格的分子筛掉,然后对前 10% 的候选分子用 32 重新精修。
Vina 还支持设置随机种子--seed。如果你希望结果可以复现,一定要加上这个参数。相同的种子、相同输入、相同参数,能得到完全一致的对接结果。这是发文章和团队协作时非常重要的一个操作,很多人会忽略。
3.3 打分函数与输出模式
Vina 默认的打分函数不需要手动指定,它会自动选择。官方文档里提到的 Vina 1.2 版本默认打分模型已经比旧版有所优化,但本质仍是基于经验的打分,不包含显式的溶剂化模型。这意味着什么?意味着 Vina 打分值不能等同于实验结合自由能。你只能说分数更负的分子比别的分子排名更靠前,不能说一个 -9.5 kcal/mol 的分子一定比 -8.5 kcal/mol 的酮强十倍。
输出模式通过--num_modes控制,默认是 9。Vina 会给出 9 个低能量构象,按打分值排序,同时输出每个构象相对于第一个构象的 RMSD。RMSD 越小,说明这两个构象越接近。如果你发现后几个模式的 RMSD 超过 2 埃,但打分只差了 0.2 kcal/mol,说明配体在这个口袋里确实存在多种相近能量构象,不能只盯着第一个模式看。
我强烈建议把 Vina 的命令日志完整保存下来。日志里有每个模式的打分值、RMSD、搜索空间大小,这些信息后期写报告或复现实验时都用得上。运行 Vina 时别忘了加--log参数:
vina --receptor receptor.pdbqt --ligand ligand.pdbqt \ --center_x 10.5 --center_y 20.3 --center_z -5.0 \ --size_x 20 --size_y 20 --size_z 20 \ --exhaustiveness 32 --num_modes 9 --seed 42 \ --out vina_out.pdbqt --log vina_log.txt上面的坐标只是一个示例,你需要替换成自己用脚本算出来的实际值。
4. 7VU6完整实战流程
4.1 受体预处理:去水、合并非极性氢、加电荷
现在正式进入 7VU6 的对接实操。假设你已经下载了 7vu6.pdb,并用文本编辑器确认里面含有共晶配体。先处理受体。
首先,打开 PyMOL,用remove solvent去掉结晶水。水分子在晶体结构里通常代表结晶溶剂,不是结合所必需。如果某个水分子恰好处在配体和蛋白之间,形成水桥相互作用,保留它会很麻烦,因为 Vina 的默认对接不支持受体水的移动。所以我的建议是:第一轮对接全部去水,如果后续发现某个水分子对结合模式影响很大,再用受体内水策略重跑。
接下来用 MGLTools 工具:
python prepare_receptor4.py -r 7vu6.pdb -A hydrogens -U nphs_lps_waters -o receptor.pdbqt这一步会完成三件事:补氢、把非极性氢合并到重原子上、计算 Gasteiger 电荷。Gasteiger 电荷是 AutoDock 系列的标准电荷模型,精度不高,但对虚拟筛选足够稳定。有些教程推荐用 AM1-BCC 等更精确的电荷,但 Vina 本身不直接用这些电荷参与打分,所以没必要过度纠结。
处理完以后,检查 receptor.pdbqt 的末尾是不是有ROOT、BRANCH这样的键信息。受体文件一般没有可旋转键,但有些非标准残基会被误判成配体片段,输出一堆 torsdof。如果出现这种情况,说明 PDB 文件里有非标准残基,你需要手动删掉或修补。
4.2 配体预处理:提取共晶配体并转PDBQT
7VU6 结构里如果有共晶配体,最好直接用这个配体来跑一次“再对接验证”(redocking)。再对接的意思是,把配体从复合物结构中提取出来,再对接到原结合位点,看能不能重现实验构象。如果能重现,说明你的对接流程和参数设置是可靠的,后续再对接新设计的分子才可信。
在 PyMOL 里提取配体很简单:
select ligand, resn LIG save ligand.pdb, sele如果配体残基名不是 LIG,先用show sequence查看链上的 HETATM,找到正确残基名再提取。
拿到 ligand.pdb 后,先加氢。如果直接从 PDB 晶体结构提取配体,通常没有氢原子,必须让工具补全。加氢时要注意质子化状态:pH 7.4 环境下,羧酸一般去质子化,氨基一般质子化。Open Babel 的通用做法是:
obabel ligand.pdb -O ligand.mol2 -p 7.4 --partialcharge gasteiger然后再转成 PDBQT:
obabel ligand.mol2 -O ligand.pdbqt --partialcharge gasteiger如果你用的是 prepare_ligand4.py,直接一步到位:
python prepare_ligand4.py -l ligand.pdb -o ligand.pdbqt这个脚本会自动判断可旋转键并生成 ROOT/BRANCH 结构。打开配体 PDBQT 文件,检查末尾的TORSDOF数值,它表示可旋转键数目。对接精度和可旋转键数目直接相关,可旋转键越多,搜索空间越复杂,结果越难稳定。如果配体有十几个可旋转键,建议把不必要的柔性链固定,或者增加 exhaustiveness。
4.3 运行Vina对接
受体和配体都准备好后,就可以跑 Vina 了。先算中心坐标,用前面那份 Python 脚本,把共晶配体坐标文件作为输入:
python calc_center.py ligand.pdb输出结果类似:
10.523 20.312 -5.078然后设置盒子尺寸。以这个配体为中心,20 埃边长通常够用。如果你想更稳一点,先跑一次盲对接,把整个蛋白包进去,看看预测的位点是否与共晶配体位置一致。如果一致,说明体系正常;如果不一致,就要查预处理哪里出了问题。
运行 Vina:
vina --receptor receptor.pdbqt --ligand ligand.pdbqt \ --center_x 10.523 --center_y 20.312 --center_z -5.078 \ --size_x 20 --size_y 20 --size_z 20 \ --exhaustiveness 32 --num_modes 9 --seed 42 \ --out vina_out.pdbqt --log vina_log.txtVina 的输出信息里会滚动显示搜索进度,完成后会在日志里显示如下表格:
mode | affinity | dist from best mode | (kcal/mol) | rmsd l.b.| rmsd u.b. -----+------------+----------+---------- 1 -8.7 0.000 0.000 2 -8.2 2.315 2.415 3 -8.0 3.102 3.205这是一份典型的模式列表。第一个模式打分最负,是 Vina 认为最优的结合构象。第二个模式与第一个的 RMSD 是 2.3 埃,打分值只差了 0.5 kcal/mol,说明存在一个能量相近但姿态明显不同的构象。如果做虚拟筛选,这两个构象都要看;如果做后续分子动力学,可以从第一个模式出发。
4.4 结果分析与可视化
Vina 的输出文件 vina_out.pdbqt 包含了所有模式,一个模式是一个 MODEL。可以用 vina_split 把它们拆成单独文件:
vina_split --input vina_out.pdbqt会生成 vina_out_ligand_1.pdbqt、vina_out_ligand_2.pdbqt 等文件。把它们导入 PyMOL,和受体、共晶配体一起显示。
在 PyMOL 里,先看共晶配体和你对接得到的第一个模式是否重叠。用命令:
align vina_out_ligand_1, ligand如果 RMSD 小于 2.0 埃,说明再对接成功,流程可靠。如果 RMSD 很大,先不要急着否定 Vina,检查两件事:一是共晶配体的周围残基有没有冲突,二是盒子中心是不是没放在口袋中心。很多时候,RMSD 偏大是因为你直接拿全蛋白对接,配体被卡在一个表面凹槽里,而不是真正的结合口袋。
接着看相互作用。Vina 不直接输出氢键列表,但你可以用 PyMOL 里的distance功能,或者在 ChimeraX 里用 “Contacts” 工具。氢键距离标准是 2.5 到 3.5 埃,疏水相互作用的距离通常在 3.5 到 4.5 埃。只看打分值不看相互作用,相当于只知相亲对象赚多少钱,不知道三观合不合。
这一步还有一个隐藏价值:如果对接构象里的配体带电基团朝向口袋外的溶剂区,或者亲水基团埋在疏水区,说明对接结果不合理。Vina 打分会因为原子重叠罚项而避免明显冲突,但它不会惩罚“亲水基团在疏水腔”这种热力学上不合理的状态。所以人工检查是必须的。
5. 常见问题与排查技巧
5.1 命令行报错定位和处理
新手最常见的状态是:Vina 装好了,文件也转换了,结果一执行就报错。别慌,绝大多数报错都能在五分钟内解决。我整理了一份速查表。
| 错误现象 | 常见原因 | 解决办法 |
|---|---|---|
| Error: could not open receptor | 文件路径写错或缺少 PDBQT 文件 | 检查文件名和路径,用ls确认 |
| Parse error on line xxx | PDBQT 文件格式损坏 | 重新运行预处理工具,不要手动乱改 |
| Unknown atom type: X | 特殊原子类型 Vina 不认识 | 用 prepare_ligand4.py 重新处理,或更换配体构象 |
| Search space is too small | 盒子尺寸小于配体大小 | 扩大size_x/y/z,保证配体可完全放入 |
| No modes satisfied | 打分都超出阈值 | 减少配体可旋转键,或增加 exhaustiveness |
| Reference ligand not found | 计算结果不包含参考结构 | 检查输出文件中是否有 MODEL 行 |
在这些报错里,最值得提的是第二种和第三种。很多初学者喜欢用文本编辑器手动修改 PDBQT,比如删掉某个原子行,结果把格式弄乱了。PDBQT 每一列都有严格对齐要求,除非你非常清楚格式规范,否则老老实实用工具生成,不要手改。
5.2 对接结果不合理?先排查这三个地方
如果 Vina 顺利跑完,但结果看起来不像实话,那问题往往在输入文件上。我每次遇到“怪结果”,都按这个顺序排查。
第一,受体结构有没有正确处理?如果 7VU6 的 PDB 文件里有多个蛋白链,比如二聚体或更复杂的组装体,但是你只用了其中一条链,可能导致结合口袋暴露在错误的表面上。很多时候配体结合需要两条链共同形成口袋,只取一条链等于把口袋撕掉一半。在用 prepare_receptor4.py 之前,先用 PyMOL 确认你要对接的生物组装单元是什么样的,必要时保留两个链。
第二,配体的质子化状态对不对?pH 会影响配体的电荷分布,比如一个羧酸在 pH 7.4 条件下应该是负电的,但如果你用的结构文件里它还是中性的 COOH,对接结果自然不对。Open Babel 的-p 7.4参数会处理 pH,但前提是输入文件里的原子连接关系正确。建议用 ChemDraw 画结构后转成 mol2,避免直接把晶体结构里的配体默认当作中性。
第三,盒子有没有包住整个配体?有些人的盒子大小刚好只够覆盖共晶配体的一部分,Vina 搜索时只要配体有一半在盒子外,就会产生不可理喻的结果。检查方法很简单:把盒子范围和对接结果一起在 PyMOL 里显示,看看配体是否完全处于盒子内部。Vina 日志里会输出搜索空间的中心和边长,你也可以根据这个信息判断。
5.3 提高重复性的两个技巧
最后分享两个我实际使用中非常受益的操作。
第一个是对接前先做一次能量最小化。这不是 Vina 能替你做的,Vina 干的是构象搜索,不是能量优化。PDB 晶体结构里的原子坐标虽然精度高,但仍可能存在很小的高能冲突。如果受体和共晶配体之间有一点点原子重叠,Vina 会为了避开冲突把构象强行扭到一边。我通常会先用 Open Babel 对配体做 MMFF94 力场的几何优化,再用 PyMOL 或 ChimeraX 的 minimize 模块对受体侧链做个几百步的能量优化。这样做能明显减少初始冲突对对接的干扰。
第二个是同一个体系至少跑三遍,每次换不同的随机种子。Vina 是随机算法,单次结果存在随机性,尤其是配体柔性较大时,两次运行可能给出不同的第一个模式。如果三遍跑下来的前三个模式都收敛到同一个构象,那么结果可信度就高。如果三遍结果完全不同,说明这个体系对搜索参数敏感,你需要提高 exhaustiveness,或者减少配体可旋转键。我自己的习惯是,正式文章里看到的“Vina 对接结果”最好都注明随机种子,否则其他人没法复现。
说到底,AutoDock Vina 只是分子对接的起点,不是终点。它能帮你快速筛掉一堆不可能结合的分子,给你一个结合模式的初判,但真正决定一个靶点能不能成药、一个分子有没有活性,还是要靠实验数据。希望这篇 7VU6 的实战流程能让你少走点弯路,跑出来的第一个结果就有参考价值。