smina-molecular-docking

smina 分子ドッキング

概要

smina は、柔軟なスコアリングと最小化に重点を置いた AutoDock Vina 1.1.2 のフォークです。SDF／MOL2／PDB リガンドを直接受け付けるため、手動で PDBQT に変換する必要はありません。参照リガンドからオートボックスを設定でき、6 種類の組み込みスコアリング関数と任意の --custom_scoring 項目が用意されており、原子ごとのスコア寄与を出力します。CLI 専用のため、Python からは subprocess 経由で実行します。

使用する場面

• 完全な探索を行わず、既存のポーズを再スコアリングまたは局所最小化する場合（--local_only、--minimize）
• ドッキングせずに単一ポーズの結合エネルギーを求める場合（--score_only）
• 標的クラスに合わせて調整したカスタムまたは経験的スコアリング関数でドッキングする場合
• リガンドごとに PDBQT へ変換せず、SDF／MOL2／複数リガンドを入力する場合
• 共結晶化した参照リガンドを囲むようにグリッドを自動設定する場合
• メディシナルケミストリーの類縁体設計に向けて、原子ごとのエネルギーを分解する場合（--atom_term_data）
• ノード間で効率よく並列化できるバッチ仮想スクリーニングを行う場合（リガンドごとに CLI プロセスを 1 つ実行）
• Vina の Python バインディング、Vinardo スコアリング、または Vina 1.2 の拡張フォースフィールドが必要な場合は、代わりに autodock-vina-docking を使用してください。結合部位が不明な場合は diffdock を使用してください。

前提条件

• smina バイナリ — conda-forge または SourceForge のビルド
• Python：rdkit、openbabel-wheel（またはシステムの openbabel）、prody、pandas、py3Dmol
• ADFR Suite：prepare_receptor 用（受容体の PDBQT のみ。リガンドは smina が処理します）
• データ：タンパク質（PDB／PDB ID）、SMILES／SDF／MOL2 形式のリガンド

インストール前に確認してください。pixi／conda 環境では、smina がすでに PATH に含まれていることがよくあります。command -v smina が成功する場合は、インストールを省略してください。pixi プロジェクト内では pixi run smina ... を実行します。

command -v smina || conda install -c conda-forge smina openbabel
pip install rdkit prody pandas py3Dmol
# ADFR Suite: https://ccsb.scripps.edu/adfr/downloads/

クイックスタート

参照リガンドからオートボックスを設定して、ドッキングを最初から最後まで実行します。

import subprocess

result = subprocess.run([
    "smina",
    "-r", "1hpv_receptor.pdbqt",
    "-l", "candidate.sdf",
    "--autobox_ligand", "1hpv_ref_ligand.pdb",
    "--autobox_add", "8",          # padding around reference (Å)
    "-o", "candidate_docked.sdf",
    "--exhaustiveness", "16",
    "--num_modes", "9",
    "--seed", "42",
], check=True, capture_output=True, text=True)

print(result.stdout.splitlines()[-15:])  # affinity table at stdout tail

ワークフロー

ステップ 1：受容体を準備する（PDBQT）

水分子とヘテロ原子を除去してから、ADFR Suite の prepare_receptor を実行します。

import subprocess, prody

pdb_id = "1HPV"
prody.fetchPDB(pdb_id, compressed=False)
protein = prody.parsePDB(f"{pdb_id}.pdb").select("protein")
prody.writePDB(f"{pdb_id}_protein.pdb", protein)

receptor_pdbqt = f"{pdb_id}_receptor.pdbqt"
subprocess.run([
    "prepare_receptor",
    "-r", f"{pdb_id}_protein.pdb",
    "-o", receptor_pdbqt,
    "-A", "hydrogens",
], check=True)
print(f"Receptor: {receptor_pdbqt} ({protein.numAtoms()} atoms)")

ステップ 2：リガンドを準備する（SDF）

smina は SDF を直接読み取ります。RDKit で 3D 座標を生成します。

from rdkit import Chem
from rdkit.Chem import AllChem

mol = Chem.MolFromSmiles("CC(C)(C)NC(=O)[C@@H]1CN(CCc2ccccc2)C[C@H]1O")
mol = Chem.AddHs(mol)
AllChem.EmbedMolecule(mol, randomSeed=42)
AllChem.MMFFOptimizeMolecule(mol)

w = Chem.SDWriter("candidate.sdf"); w.write(mol); w.close()
print(f"Ligand SDF: candidate.sdf ({mol.GetNumAtoms()} atoms)")

ステップ 3：オートボックス用の参照リガンドを抽出する

--autobox_ligand を使うと、参照構造からグリッドが設定されます。

import prody

ref = prody.parsePDB(f"{pdb_id}.pdb").select("hetero and not water and not ion")
if ref is None:
    raise RuntimeError("No reference ligand — supply explicit --center_x/--size_x")
prody.writePDB(f"{pdb_id}_ref_ligand.pdb", ref)
print(f"Ref ligand: {ref.numAtoms()} atoms, center {ref.getCoords().mean(axis=0).round(2)}")

ステップ 4：オートボックスでドッキングを実行する

親和性の表は標準出力に表示されるため、取得してください。

import subprocess

proc = subprocess.run([
    "smina",
    "-r", receptor_pdbqt,
    "-l", "candidate.sdf",
    "--autobox_ligand", f"{pdb_id}_ref_ligand.pdb",
    "--autobox_add", "8",
    "-o", "candidate_docked.sdf",
    "--exhaustiveness", "16",
    "--num_modes", "9",
    "--energy_range", "3",
    "--cpu", "4",
    "--seed", "42",
], check=True, capture_output=True, text=True)

for line in proc.stdout.splitlines()[-15:]:
    print(line)

ステップ 5：ポーズと親和性を解析する

親和性は SDF の <minimizedAffinity> プロパティに格納されます。

from rdkit import Chem
import pandas as pd

rows = []
for i, mol in enumerate(Chem.SDMolSupplier("candidate_docked.sdf", removeHs=False)):
    if mol is None:
        continue
    aff = float(mol.GetProp("minimizedAffinity")) if mol.HasProp("minimizedAffinity") else None
    rmsd = float(mol.GetProp("minimizedRMSD")) if mol.HasProp("minimizedRMSD") else None
    rows.append({"pose": i + 1, "affinity_kcal_mol": aff, "rmsd_to_best": rmsd})

df = pd.DataFrame(rows).sort_values("affinity_kcal_mol")
print(df.to_string(index=False))
print(f"Best: {df.iloc[0]['affinity_kcal_mol']:.2f} kcal/mol")

ステップ 6：既存のポーズを局所最小化する

--local_only を指定すると、グローバル探索を行わずに入力ポーズを精密化します。

import subprocess

subprocess.run([
    "smina",
    "-r", receptor_pdbqt,
    "-l", "candidate_pose.sdf",
    "--autobox_ligand", "candidate_pose.sdf",   # box around the pose itself
    "--autobox_add", "4",
    "-o", "candidate_min.sdf",
    "--local_only",
    "--minimize_iters", "1000",
], check=True, capture_output=True, text=True)
print("Local minimization complete: candidate_min.sdf")

ステップ 7：ドッキング複合体を可視化する

import py3Dmol

with open(f"{pdb_id}_protein.pdb") as f: rec = f.read()
with open("candidate_docked.sdf") as f: lig = f.read()

view = py3Dmol.view(width=800, height=600)
view.addModel(rec, "pdb"); view.setStyle({"model": 0}, {"cartoon": {}})
view.addModel(lig.split("$$$$")[0] + "$$$$", "sdf")
view.setStyle({"model": 1}, {"stick": {}})
view.zoomTo({"model": 1}); view.show()

ステップ 8：バッチ仮想スクリーニング

リガンドごとに smina プロセスを 1 つ実行すると、コア間またはクラスターノード間で効率よく並列化できます。

import subprocess, pandas as pd
from rdkit import Chem
from rdkit.Chem import AllChem

library = pd.DataFrame({
    "name":   ["cpd_001", "cpd_002", "cpd_003"],
    "smiles": ["CC(=O)Oc1ccccc1C(=O)O",
               "CC(C)Cc1ccc(cc1)C(C)C(=O)O",
               "OC(=O)c1ccccc1O"],
})

results = []
for _, row in library.iterrows():
    lig_sdf, out_sdf = f"{row['name']}.sdf", f"{row['name']}_docked.sdf"
    mol = Chem.MolFromSmiles(row["smiles"]); mol = Chem.AddHs(mol)
    AllChem.EmbedMolecule(mol, randomSeed=42); AllChem.MMFFOptimizeMolecule(mol)
    w = Chem.SDWriter(lig_sdf); w.write(mol); w.close()

    proc = subprocess.run([
        "smina", "-r", receptor_pdbqt, "-l", lig_sdf,
        "--autobox_ligand", f"{pdb_id}_ref_ligand.pdb", "--autobox_add", "8",
        "-o", out_sdf, "--exhaustiveness", "8", "--num_modes", "1",
        "--cpu", "2", "--seed", "42",
    ], capture_output=True, text=True)

    if proc.returncode != 0:
        results.append({"name": row["name"], "affinity_kcal_mol": None})
        continue
    docked = next(iter(Chem.SDMolSupplier(out_sdf, removeHs=False)))
    aff = float(docked.GetProp("minimizedAffinity")) if docked and docked.HasProp("minimizedAffinity") else None
    results.append({"name": row["name"], "affinity_kcal_mol": aff})

ranked = pd.DataFrame(results).sort_values("affinity_kcal_mol")
ranked.to_csv("screening_results.csv", index=False)
print(ranked.to_string(index=False))

主なパラメーター

パラメーター	デフォルト	範囲／オプション	効果
--exhaustiveness	8	1-128	MC 探索の試行量。実運用では 16-32、論文発表向けには 64 以上
--num_modes	9	1-20	出力 SDF に書き込むポーズ数
--energy_range	3	1-5 kcal/mol	出力する最良ポーズからの最大 ΔE
--scoring	default	default、vina、vinardo、dkoes_scoring、dkoes_scoring_old、ad4_scoring	組み込みスコアリング（default は約 Vina 相当）
--custom_scoring	—	path	ユーザー定義の重み付き項目。--scoring を上書きします
--autobox_ligand	—	PDB/SDF path	参照構造からグリッドを設定
--autobox_add	4	2-12 Å	オートボックス範囲の周囲に追加するパディング
--center_x/y/z、--size_x/y/z	—	Å	手動で設定するグリッド（オートボックスの代替）
--local_only	off	flag	グローバル探索をスキップし、局所最適化のみ実行
--score_only	off	flag	最小化や探索を行わずにエネルギーを算出
--minimize	off	flag	スコアリング探索を行わずにエネルギーを最小化
--minimize_iters	0 (auto)	0-100000	局所最小化の反復回数
--atom_term_data	off	flag	出力に原子ごとのスコア寄与を含める
--cpu	all	1-N	ジョブごとのスレッド数
--seed	random	int	再現性を確保するための RNG シード
--no_lig	off	flag	ベースラインとして受容体のみをスコアリング

よく使うレシピ

レシピ：スコアのみ（単一点エネルギー）

再ドッキングせずに、スコアリング関数間でポーズのエネルギーを比較します。

import subprocess

proc = subprocess.run([
    "smina", "-r", "receptor.pdbqt", "-l", "pose.sdf",
    "--score_only", "--scoring", "vinardo",   # try vina, vinardo, dkoes_scoring
], capture_output=True, text=True, check=True)

for line in proc.stdout.splitlines():
    if line.startswith("Affinity:"):
        print(line)

レシピ：カスタムスコアリング関数

ターゲットクラスに合わせて調整した経験的スコアリングです。dkoes の項、またはカスタムフィットした項のセットを再現します。

import subprocess

with open("my_scoring.txt", "w") as f:
    f.write("""\
-0.035579    gauss(o=0,_w=0.5,_c=8)
-0.005156    gauss(o=3,_w=2,_c=8)
0.840245     repulsion(o=0,_c=8)
-0.035069    hydrophobic(g=0.5,_b=1.5,_c=8)
-0.587439    non_dir_h_bond(g=-0.7,_b=0,_c=8)
1.923        num_tors_div
""")

subprocess.run([
    "smina", "-r", "receptor.pdbqt", "-l", "candidate.sdf",
    "--custom_scoring", "my_scoring.txt",
    "--autobox_ligand", "ref.pdb", "--autobox_add", "8",
    "-o", "custom_docked.sdf", "--exhaustiveness", "16",
], check=True)

レシピ：原子ごとのスコア分解

どのリガンド原子が結合に寄与しているかを特定します。

import subprocess
from rdkit import Chem

subprocess.run([
    "smina", "-r", "receptor.pdbqt", "-l", "pose.sdf",
    "--score_only", "--atom_term_data", "-o", "pose_decomp.sdf",
], check=True)

mol = next(iter(Chem.SDMolSupplier("pose_decomp.sdf", removeHs=False)))
for prop in mol.GetPropNames():
    if "atom_term" in prop:
        print(f"{prop}: {mol.GetProp(prop)[:120]}")

レシピ：再ドッキングによる検証

共結晶化したリガンドを再ドッキングし、RMSD < 2.0 Å を確認します。

import subprocess
from rdkit import Chem
from rdkit.Chem import AllChem

subprocess.run([
    "smina", "-r", "receptor.pdbqt", "-l", "ref_ligand.sdf",
    "--autobox_ligand", "ref_ligand.sdf", "--autobox_add", "8",
    "-o", "redocked.sdf", "--exhaustiveness", "32",
    "--num_modes", "1", "--seed", "42",
], check=True)

ref = next(iter(Chem.SDMolSupplier("ref_ligand.sdf", removeHs=False)))
docked = next(iter(Chem.SDMolSupplier("redocked.sdf", removeHs=False)))
rmsd = AllChem.GetBestRMS(ref, docked)
print(f"Re-docking RMSD: {rmsd:.2f} Å  ->  {'PASS' if rmsd < 2.0 else 'FAIL'}")

期待される出力

• *_docked.sdf — マルチモデル SDF。ポーズごとに1分子を格納し、結合親和性は <minimizedAffinity>、最良ポーズに対する RMSD は <minimizedRMSD> に格納
• screening_results.csv — 名前、SMILES、結合親和性（kcal/mol）
• *_receptor.pdbqt — 調製済み受容体
• stdout — 結合親和性テーブル：ポーズ、結合親和性（kcal/mol）、RMSD の下限／上限

トラブルシューティング

問題	原因	解決策
smina: command not found	バイナリが PATH にない	conda install -c conda-forge smina を実行します。pixi 環境では pixi run smina を使用します
Parse error on line ...（リガンド）	SDF／MOL2 の形式が不正	AddHs と EmbedMolecule の実行後、RDKit の SDWriter で再エクスポートします
Could not figure out box dimensions	--autobox_ligand と --center_x/--size_x の両方が指定されていない	オートボックスの参照または明示的なグリッドを指定します
結合親和性が非常に正の値（>0）	リガンドがグリッドの外にある、または入力に立体障害がある	--autobox_add の値を上げ、参照ポーズが結合部位内にあることを確認します
--local_only が入力を変更せずに返す	一部のビルドでは --minimize_iters 0 は「最小化しない」という意味です	--minimize_iters 1000 を明示的に設定します
実行ごとに同一のポーズになる	--seed が固定されており、--exhaustiveness の値が低い	32 以上に上げ、アンサンブルには --seed を変えます
受容体の PDBQT 生成に失敗する	prepare_receptor が PATH にない	ADFR Suite をインストールし、<adfr>/bin を PATH に追加します
ドッキング後の SDF が空	smina が警告なく中断した	capture_output を削除して proc.stderr を確認し、受容体の電荷をチェックします
--score_only と --local_only でスコアが異なる	ローカル最適化でポーズが移動する	想定どおりです。固定ポーズには score-only、精密化したエネルギーには local-only を使用します
--custom_scoring ファイルが拒否される	項の名前の誤り、または重みの列がない	各行は <weight><whitespace><term_name(args)> の形式です。空白の形式は厳密です
原子ごとの分解結果が空	--score_only ではなく --minimize／ドッキングと併用した	--atom_term_data は --score_only との併用が最も確実です

参考資料

• SourceForge の smina — オリジナルの配布元
• smina の GitHub ミラー — メンテナンスされているビルド
• Koes、Baumgartner、Camacho（2013） — smina と dkoes スコアリングに関する論文、J Chem Inf Model
• Trott と Olson（2010） — 上流の Vina 論文、J Comput Chem
• ADFR Suite — 受容体 PDBQT 用の prepare_receptor
• RDKit のドキュメント — リガンドの3D生成と SDF 入出力
• Open Babel — リガンド形式変換の代替手段