autodock-vina-docking

AutoDock Vina 分子对接

概述

AutoDock Vina 是速度最快、应用最广泛的开源分子对接引擎之一，可用于预测蛋白质–配体的结合模式和亲和力。本技能涵盖完整的 Python 流程：从 PDB 准备受体，使用 Meeko 和 RDKit 从 SMILES/SDF 准备配体，定义搜索框，执行对接，分析构象，以及通过批量虚拟筛选识别命中化合物。

适用场景

• 预测小分子与蛋白质靶标的结合构象
• 估算相对结合亲和力（kcal/mol），以对配体进行排序
• 针对靶标受体筛选化合物库
• 通过重新对接共结晶配体验证对接方案
• 无需中间文件，直接从 SMILES 字符串准备对接输入
• 在构效关系研究中比较类似物的结合模式
• 为分子动力学模拟生成初始构象
• 如果结合位点未知，盲对接请改用 DiffDock；也可使用 GNINA 作为采用 CNN 评分的替代方案

前置条件

• Python 包：vina、meeko、rdkit、prody（用于获取 PDB）、py3Dmol（用于可视化）
• 外部工具：ADFR Suite（提供用于 PDBQT 转换的 prepare_receptor）——从 https://ccsb.scripps.edu/adfr/downloads/ 下载
• 数据要求：蛋白质结构（PDB 文件或 PDB ID），配体可采用 SMILES、SDF 或 MOL2 格式
• 环境：建议使用 Python 3.8+、Linux 或 macOS

pip install vina meeko rdkit-pypi prody py3Dmol
# ADFR Suite must be installed separately for prepare_receptor

工作流

步骤 1：获取并准备受体

下载蛋白质结构，移除水分子和杂原子，并转换为 PDBQT 格式。

import prody
import subprocess
from pathlib import Path

# Download PDB structure (example: HIV-1 protease, PDB 1HPV)
pdb_id = "1HPV"
pdb_file = f"{pdb_id}.pdb"
prody.fetchPDB(pdb_id, compressed=False)

# Extract protein chain only (remove water and ligands)
structure = prody.parsePDB(pdb_file)
protein = structure.select("protein")
prody.writePDB(f"{pdb_id}_protein.pdb", protein)
print(f"Protein atoms: {protein.numAtoms()}")

# Convert to PDBQT using ADFR Suite's prepare_receptor
receptor_pdbqt = f"{pdb_id}_receptor.pdbqt"
subprocess.run([
    "prepare_receptor",
    "-r", f"{pdb_id}_protein.pdb",
    "-o", receptor_pdbqt,
    "-A", "hydrogens",  # add hydrogens
], check=True)
print(f"Receptor PDBQT: {receptor_pdbqt}")

步骤 2：确定结合位点

以已知结合位点或共结晶配体为中心，定义对接搜索框。

import numpy as np

# Option A: Center on co-crystallized ligand coordinates
structure = prody.parsePDB(pdb_file)
ligand = structure.select("hetero and not water and not ion")
if ligand is not None:
    center = ligand.getCoords().mean(axis=0)
    # Box size = ligand extent + padding
    extent = ligand.getCoords().max(axis=0) - ligand.getCoords().min(axis=0)
    box_size = extent + 10.0  # 10 Å padding on each side
    print(f"Binding site center: {center}")
    print(f"Box size: {box_size}")
else:
    # Option B: Manual coordinates (from literature or visual inspection)
    center = np.array([15.0, 54.0, 17.0])
    box_size = np.array([25.0, 25.0, 25.0])
    print(f"Using manual box: center={center}, size={box_size}")

步骤 3：从 SMILES 准备配体

使用 RDKit 生成 3D 坐标，并使用 Meeko 转换为 PDBQT 格式。

from rdkit import Chem
from rdkit.Chem import AllChem
from meeko import MoleculePreparation, PDBQTWriterLegacy

# Define ligand (example: Indinavir, HIV protease inhibitor)
smiles = "CC(C)(C)NC(=O)[C@@H]1CN(CCc2ccccc2)C[C@H]1O"
mol_name = "indinavir_analog"

# Generate 3D coordinates with RDKit
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
AllChem.EmbedMolecule(mol, randomSeed=42)
AllChem.MMFFOptimizeMolecule(mol)

# Convert to PDBQT using Meeko
preparator = MoleculePreparation()
mol_setups = preparator.prepare(mol)

# Write PDBQT string (for Vina API) or file
pdbqt_string = PDBQTWriterLegacy.write_string(mol_setups[0])[0]
ligand_pdbqt = f"{mol_name}.pdbqt"
with open(ligand_pdbqt, "w") as f:
    f.write(pdbqt_string)
print(f"Ligand PDBQT: {ligand_pdbqt} ({mol.GetNumAtoms()} atoms)")

步骤 4：运行对接

初始化 AutoDock Vina，配置搜索空间并执行对接。

from vina import Vina

# Initialize Vina
v = Vina(sf_name="vina", cpu=4)

# Load receptor and ligand
v.set_receptor(receptor_pdbqt)
v.set_ligand_from_file(ligand_pdbqt)

# Define search space
v.compute_vina_maps(
    center=center.tolist(),
    box_size=box_size.tolist(),
)

# Run docking
v.dock(
    exhaustiveness=32,   # Higher = more thorough (default 8)
    n_poses=10,          # Number of output poses
)

# Write output poses
output_file = f"{mol_name}_docked.pdbqt"
v.write_poses(output_file, n_poses=10, overwrite=True)
print(f"Docking complete. Poses written to {output_file}")

步骤 5：分析对接结果

使用 Vina 内置 API 从对接构象中提取结合能和 RMSD 值。

# Method A: Vina built-in energies (most reliable)
energies = v.energies(n_poses=10)
# Each row: [total, inter, intra, torsions, intra_best_pose]
print(f"{'Pose':<6} {'Total (kcal/mol)':<18} {'Inter':<10} {'Intra':<10}")
print("-" * 44)
for i, e in enumerate(energies):
    print(f"{i+1:<6} {e[0]:<18.2f} {e[1]:<10.2f} {e[2]:<10.2f}")

print(f"\nBest pose: {energies[0][0]:.2f} kcal/mol")

# Method B: Parse PDBQT output with Meeko (for RDKit conversion)
from meeko import PDBQTMolecule, RDKitMolCreate
pdbqt_mol = PDBQTMolecule.from_file(output_file)
best_pose_rdkit = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)[0]
print(f"Best pose converted to RDKit mol: {best_pose_rdkit.GetNumAtoms()} atoms")

步骤 6：可视化对接结果

使用 py3Dmol 生成对接复合物的 3D 可视化。

import py3Dmol

# Load receptor and best docked pose
with open(f"{pdb_id}_protein.pdb") as f:
    receptor_pdb = f.read()
with open(output_file) as f:
    docked_pdbqt = f.read()

# Create 3D viewer
view = py3Dmol.view(width=800, height=600)
view.addModel(receptor_pdb, "pdb")
view.setStyle({"model": 0}, {"cartoon": {"color": "lightgrey"}})

# Add docked ligand (first model only)
first_model = docked_pdbqt.split("ENDMDL")[0] + "ENDMDL"
view.addModel(first_model, "pdb")
view.setStyle({"model": 1}, {"stick": {"colorscheme": "greenCarbon"}})

# Zoom to ligand
view.zoomTo({"model": 1})
view.show()
# In Jupyter: displays interactive 3D view
# To save: view.png() or view.write_html("docking_result.html")
print("3D visualization rendered")

步骤 7：批量虚拟筛选

使用同一受体筛选化合物库。

import pandas as pd

# Define compound library
compounds = pd.DataFrame({
    "name": ["cpd_001", "cpd_002", "cpd_003", "cpd_004", "cpd_005"],
    "smiles": [
        "CC(=O)Oc1ccccc1C(=O)O",          # Aspirin
        "CC(C)Cc1ccc(cc1)C(C)C(=O)O",      # Ibuprofen
        "OC(=O)c1ccccc1O",                  # Salicylic acid
        "CC12CCC3C(CCC4CC(=O)CCC34C)C1CCC2O",  # Testosterone
        "c1ccc2c(c1)cc1ccc3cccc4ccc2c1c34",     # Pyrene
    ],
})

# Screen all compounds
results = []
for _, row in compounds.iterrows():
    try:
        # Prepare ligand
        mol = Chem.MolFromSmiles(row["smiles"])
        mol = Chem.AddHs(mol)
        AllChem.EmbedMolecule(mol, randomSeed=42)
        AllChem.MMFFOptimizeMolecule(mol)

        mol_setups = preparator.prepare(mol)
        pdbqt_str = PDBQTWriterLegacy.write_string(mol_setups[0])[0]

        # Dock
        v_screen = Vina(sf_name="vina", cpu=2)
        v_screen.set_receptor(receptor_pdbqt)
        v_screen.set_ligand_from_string(pdbqt_str)
        v_screen.compute_vina_maps(center=center.tolist(), box_size=box_size.tolist())
        v_screen.dock(exhaustiveness=16, n_poses=1)

        energy = v_screen.energies(n_poses=1)[0][0]
        results.append({"name": row["name"], "smiles": row["smiles"], "energy_kcal": energy})
        print(f"  {row['name']}: {energy:.2f} kcal/mol")

    except Exception as e:
        results.append({"name": row["name"], "smiles": row["smiles"], "energy_kcal": None})
        print(f"  {row['name']}: FAILED ({e})")

# Rank by binding energy
results_df = pd.DataFrame(results).sort_values("energy_kcal")
results_df.to_csv("screening_results.csv", index=False)
print(f"\nTop hits:\n{results_df.head()}")

步骤 8：保存并导出

import os
os.makedirs("results", exist_ok=True)

# Save summary
results_df.to_csv("results/screening_results.csv", index=False)

# Save best poses for top hits
for _, row in results_df.head(3).iterrows():
    print(f"Top hit: {row['name']} → {row['energy_kcal']:.2f} kcal/mol")

print("Virtual screening complete. Results in results/screening_results.csv")

关键参数

参数	默认值	范围 / 选项	作用
exhaustiveness	8	8-128	搜索彻底程度；发表论文时建议设为 32+
n_poses	9	1-20	输出的结合构象数量
energy_range	3.0	1.0-5.0	纳入的构象与最佳构象之间的最大能量差（kcal/mol）
sf_name	"vina"	"vina"、"ad4"、"vinardo"	选择评分函数
cpu	all	1-N	用于对接的 CPU 核心数
box_size (xyz)	—	每边 15-30 Å	搜索空间尺寸；必须涵盖结合位点，并留出 5-10Å 的边缘空间
center (xyz)	—	结合位点坐标	搜索框中心
randomSeed (RDKit)	random	任意 int	生成可复现的 3D 构象
padding (box)	10.0 Å	5.0-15.0 Å	定义搜索框时，在已知配体周围增加的空间

常用方案

方案：重新对接验证（同源配体对接）

适用场景：通过重新对接共结晶配体并检查 RMSD < 2.0 Å 来验证方案。

from rdkit.Chem import AllChem, rdMolAlign

# Extract co-crystallized ligand from PDB
ref_ligand = Chem.MolFromPDBFile(f"{pdb_id}_ligand.pdb", removeHs=False)

# Dock the same ligand
# ... (use steps 3-4 above with the extracted ligand)

# Calculate RMSD between docked pose and crystal structure
rmsd = AllChem.GetBestRMS(ref_ligand, best_pose_rdkit)
print(f"Re-docking RMSD: {rmsd:.2f} Å")
print(f"Validation: {'PASS' if rmsd < 2.0 else 'FAIL'} (threshold: 2.0 Å)")

方案：柔性受体对接

适用场景：结合位点的关键残基需要具有构象自由度（例如，诱导契合）。

# Prepare receptor with flexible sidechains (using ADFR Suite)
# prepare_receptor -r protein.pdb -o rigid.pdbqt -A hydrogens
# prepare_flexreceptor -r rigid.pdbqt -s "A:ARG8,A:ASP25,A:ILE50"

v_flex = Vina(sf_name="vina", cpu=4)
v_flex.set_receptor("rigid.pdbqt", "flex.pdbqt")  # rigid + flexible parts
v_flex.set_ligand_from_file(ligand_pdbqt)
v_flex.compute_vina_maps(center=center.tolist(), box_size=box_size.tolist())
v_flex.dock(exhaustiveness=64, n_poses=10)
v_flex.write_poses("docked_flex.pdbqt", n_poses=5, overwrite=True)

方案：仅评分（不进行对接）

适用场景：无需运行完整搜索，即可评估已预先定位的配体的结合能。

v_score = Vina(sf_name="vina")
v_score.set_receptor(receptor_pdbqt)
v_score.set_ligand_from_file("pre_positioned_ligand.pdbqt")
v_score.compute_vina_maps(center=center.tolist(), box_size=box_size.tolist())

# Score current pose
energy = v_score.score()
print(f"Score: {energy[0]:.2f} kcal/mol")

# Local minimization
energy_min = v_score.optimize()
print(f"After local optimization: {energy_min[0]:.2f} kcal/mol")
v_score.write_pose("minimized.pdbqt", overwrite=True)

方案：比较多种评分函数

适用场景：通过共识评分提高对接结果的可信度。

scoring_results = {}
for sf in ["vina", "vinardo", "ad4"]:
    v_sf = Vina(sf_name=sf, cpu=2)
    v_sf.set_receptor(receptor_pdbqt)
    v_sf.set_ligand_from_file(ligand_pdbqt)
    v_sf.compute_vina_maps(center=center.tolist(), box_size=box_size.tolist())
    v_sf.dock(exhaustiveness=16, n_poses=1)
    scoring_results[sf] = v_sf.energies(n_poses=1)[0][0]

for sf, energy in scoring_results.items():
    print(f"  {sf}: {energy:.2f} kcal/mol")

预期输出

• {name}_docked.pdbqt — PDBQT 格式的对接构象，头部包含结合能
• results/screening_results.csv — 虚拟筛选结果：化合物名称、SMILES、结合能（kcal/mol）
• {pdb_id}_receptor.pdbqt — PDBQT 格式的受体准备文件
• 图：3D 对接可视化（交互式 py3Dmol 或静态图像）
• 控制台输出：按排名列出的结合能和每个构象的 RMSD 值

故障排除

问题	原因	解决方案
找不到 prepare_receptor	ADFR Suite 不在 PATH 中	将 ADFR Suite bin 添加到 $PATH，或使用完整路径
RuntimeError: receptor not set	忘记调用 set_receptor	对接前调用 v.set_receptor(pdbqt_file)
对接分数非常高（>0）	配体在搜索框外或几何结构不佳	检查搜索框中心/尺寸是否覆盖结合位点；确认 3D 坐标
所有构象都相同	exhaustiveness 过低	增至 32-64，以确保可靠采样
Meeko MoleculePreparation 错误	输入分子缺少氢原子	使用 Meeko 前务必调用 Chem.AddHs(mol)
重对接时 RMSD > 2Å	搜索框过小或中心位置错误	将搜索框扩大 5Å；确认中心位于共结晶配体处
EmbedMolecule 返回 -1	RDKit 未能生成 3D 坐标	使用 AllChem.EmbedMolecule(mol, maxAttempts=1000)，或尝试 useRandomCoords=True
筛选速度慢（每个化合物 >1 分钟）	exhaustiveness 值高且搜索框较大	筛选时将 exhaustiveness 降至 8-16；缩小搜索框
找不到 PDBQTWriterLegacy	Meeko 版本较旧	pip install meeko>=0.5 — API 已从 write_pdbqt_string 更改
多次运行的能量不一致	搜索过程不确定	在 v.dock(seed=42) 中设置 seed 参数，以确保结果可复现

附带资源

此技能包含用于深入查阅的参考文件。需要时再阅读。

references/receptor_preparation_guide.md

受体准备详细指南：处理缺失残基、质子化状态（依赖 pH）、金属离子、辅因子和多链复合物。提供决策树，说明何时使用 PDB2PQR、PROPKA 或手动质子化。

references/scoring_functions_comparison.md

Vina、Vinardo 和 AD4 评分函数的对比：准确度基准、速度权衡，以及针对不同靶标类别（激酶、蛋白酶、GPCR、核受体）的建议。

参考资料

• AutoDock Vina 文档 — 官方文档和 Python API（Apache-2.0）
• Eberhardt et al. (2021) — “AutoDock Vina 1.2.0: New Docking Methods, Expanded Force Field, and Python Bindings”，J Chem Inf Model
• Meeko：为 AutoDock 准备小分子 — 使用 RDKit 准备配体 PDBQT 文件（LGPL）
• ADFR Suite — Scripps 提供的受体准备工具（免费学术许可）
• Forli et al. (2016) — “Computational protein-ligand docking and virtual drug screening with the AutoDock suite”，Nat Protoc
• RDKit 文档 — 用于配体处理的化学信息学工具包（BSD）