autodock-vina-docking

AutoDock Vina 分子ドッキング

概要

AutoDock Vina は、タンパク質とリガンドの結合様式および親和性を予測する、最速かつ最も広く利用されているオープンソースの分子ドッキングエンジンの一つです。このスキルでは、PDB からの受容体調製、Meeko と RDKit による SMILES／SDF からのリガンド調製、探索ボックスの定義、ドッキングの実行、ポーズ解析、ヒット化合物の特定に向けたバッチバーチャルスクリーニングまで、Python ベースの一連のパイプラインを解説します。

使用する場面

• タンパク質標的に対する小分子の結合ポーズを予測する場合
• リガンドの順位付けのため、相対結合親和性（kcal/mol）を推定する場合
• 標的受容体に対して化合物ライブラリのバーチャルスクリーニングを行う場合
• 共結晶化リガンドを再ドッキングして、ドッキングプロトコルを検証する場合
• 中間ファイルを作成せずに、SMILES 文字列からドッキング用入力データを準備する場合
• 構造活性相関の研究で類縁体の結合様式を比較する場合
• 分子動力学シミュレーション用の初期ポーズを生成する場合
• 結合部位が不明な場合は、ブラインドドッキングには DiffDock を代わりに使用してください。CNN スコアリングを使う代替手段として GNINA も使用できます

前提条件

• 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)	—	1 辺あたり 15-30 Å	探索空間の寸法。結合部位に 5-10Å の余裕を加えた範囲を含める必要があります
center (xyz)	—	結合部位の座標	探索ボックスの中心
randomSeed (RDKit)	random	任意の整数	再現可能な 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 エラー	入力 mol に水素がない	Meeko の前に必ず Chem.AddHs(mol) を呼び出します
再ドッキングで RMSD > 2Å	ボックスが小さすぎる、または中心が誤っている	ボックスを 5Å 拡大し、共結晶リガンド上に中心があることを確認します
EmbedMolecule が -1 を返す	RDKit による 3D 座標の生成に失敗	AllChem.EmbedMolecule(mol, maxAttempts=1000) を使用するか、useRandomCoords=True を試します
スクリーニングが遅い（1 化合物あたり >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 documentation — 公式ドキュメントと 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: Preparation of small molecules for 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 documentation — リガンド処理用のケモインフォマティクスツールキット（BSD）