本記事はMacでのin silico screeningの方法を書いた記事です。AutoDock Vinaによるin silico screeningはMacで行うと、環境構築が難しいですが、AutoDockの改良版であるsminaを使うと簡単にすることができます。Macの方で、in silico screeningを行いたかった方、ぜひトライしてみてください。
Mac M1, Sequoia 15.3, メモリ 16GB
自宅でできるin silico創薬の技術書を販売中
新薬探索を試したい方必読！
ITエンジニアである著者の視点から、wetな研究者からもdryの創薬研究をわかりやすく身近に感じられるように解説しています
目次
Smina とは？
インシリコスクリーニングは、数千～数百万の化合物ライブラリを仮想的に試し、ターゲットタンパク質に最も強く結合する候補を探す手法です。代表的なツールとして AutoDock Vina がありますが、 Smina はその改良版です。
AutoDock Vina は 高速かつ使いやすい ドッキングツールで、スコアリング関数は固定ですが、並列処理により効率的なスクリーニングが可能です。一方、Smina は スコアリング関数を変更できる ため、AutoDock4（AD4）やカスタム関数を用いた評価が可能です。また、Smina では --score_only 機能により、ドッキングせずに 結合エネルギーを計算するだけ の処理ができ、大規模スクリーニングの高速化に適しています。さらに、--autobox_ligand を利用すると、リガンドサイズに基づいて自動的にグリッドボックスを設定でき、スクリーニングの効率が向上します。
またAutoDock vinaで必要なpdbqt変換がなしで、済むので、準備が楽です。
それでは実際にSminaを試していきたいと思います！
タンパク質準備
UCFS chimeraで行います。最近ではChimera Xというバージョンが上がっていますが、こちらだとエネルギー最小化ができないそうなので、古いバージョンで行います。
こちらについては過去の記事を参考にしてください。この記事ではタンパク質は糖尿病の標的タンパク質DPPIV（PDBID: 2I78)を使います。
インシリコ創薬の分子ドッキングで必要なタンパク質の準備を論文に沿ってやってみる。【in silico創薬】
拙著
エンジニアとWET研究者のためのコンピュータ創薬入門 ~コンピュータ創薬を体験して、新しい技術を習得しよう~
自宅でできるin silico創薬【新薬探索編】 ~初心者でもできる”おうち創薬”~
では本記事と同タンパク質でタンパク質の準備をしております。よろしければご参照ください。
リガンドの準備
今回はすぐに実験に使えるライブラリとして、selleckchem bioactive library Iを用います。以下を参照してください。
化合物ライブラリまとめ【in silico創薬】【in silicoスクリーニング】 – LabCode
上記から入手したsdfファイルをrule of fiveやPAINSフィルターでフィルタリングします。興味があれば、以下の記事をご参照ください。やらなくても大丈夫です。
【化合物スクリーニング】RDkitを使った化合物のフィルタリング【in silico創薬】 – LabCode
化合物の分割
SDFファイルに対して、まずターミナルで以下を実行していきます。
obabel filtered_compounds.sdf -O drug.sdf -m
このコマンドは、Open Babel (obabel) を使って filtered_compounds.sdf ファイルを drug.sdf に変換し、個々の分子を分割する ためのものです。
filtered_compounds.sdfの分子データを読み込む.sdf(Structure Data File) は、化合物の3D構造や付随情報を格納するファイル形式。
drug.sdfに変換・出力する- 変換結果を
drug.sdfという名前で保存する。
- 変換結果を
mオプションで分子ごとに個別ファイルを作成- 通常、
.sdf形式は 複数の分子を1つのファイルにまとめる 仕様。 mを指定すると、各分子が**drug1.sdf,drug2.sdf,drug3.sdf, …** のように 個別のファイル として出力される。
- 通常、
続いて、Pythonで以下のコードを実行していってください。環境構築はこちらを参考にしてみてください。
【化合物スクリーニング】RDkitを使った化合物のフィルタリング【in silico創薬】 – LabCode
水分子、イオンの除去
リガンドに一緒にある水分子やイオンを除いています。
import glob
from rdkit import Chem
from rdkit.Chem import SDMolSupplier, SDWriter
# 除去する原子と水分子
REMOVE_ATOMS = {"Cl", "Na", "K", "Ca", "Mg"}
REMOVE_SMARTS = ["[OH2]"] # H2O の SMARTS パターン
def remove_unwanted_atoms(mol):
"""Cl, Na, K, Ca, Mg の単独原子を削除"""
if mol is None:
return None
edit_mol = Chem.RWMol(mol)
remove_idxs = [atom.GetIdx() for atom in mol.GetAtoms() if atom.GetSymbol() in REMOVE_ATOMS]
for idx in sorted(remove_idxs, reverse=True):
edit_mol.RemoveAtom(idx)
return edit_mol.GetMol() # RWMol → ROMol に変換
def remove_water_molecules(mol):
"""H2O の SMARTS パターンに一致する分子を削除"""
if mol is None:
return None
for smarts in REMOVE_SMARTS:
pattern = Chem.MolFromSmarts(smarts)
if mol.HasSubstructMatch(pattern):
return None # 水分子なら削除
return mol
def process_sdf(input_sdf, output_sdf):
"""SDF ファイルを処理し、不要な原子や水分子を削除"""
supplier = SDMolSupplier(input_sdf, removeHs=False)
writer = SDWriter(output_sdf)
for mol in supplier:
if mol is not None:
mol = remove_unwanted_atoms(mol) # Cl, Na, K など削除
mol = remove_water_molecules(mol) # 水分子削除
if mol is not None:
try:
# **Kekulé 化のエラーを防ぐ**
Chem.SetAromaticity(mol, Chem.AromaticityModel.AROMATICITY_MDL)
Chem.SanitizeMol(mol, Chem.SanitizeFlags.SANITIZE_KEKULIZE)
# **RWMol → ROMol に変換後、Kekulé 化をスキップして保存**
writer.write(mol)
except Chem.KekulizeException:
print(f"Warning: Failed to kekulize {input_sdf}. Skipping this molecule.")
writer.close()
print(f"Processed: {output_sdf}")
# カレントディレクトリ内の全 .sdf ファイルを処理
input_sdf_files = glob.glob("*.sdf") # フォルダ内の .sdf ファイルを取得
output_sdf_files = [f"clean_{fname}" for fname in input_sdf_files]
for input_sdf, output_sdf in zip(input_sdf_files, output_sdf_files):
process_sdf(input_sdf, output_sdf)
このコードは、RDKit を用いて SDF ファイル（化学構造データ）を処理し、不要な原子や水分子を削除する スクリプトです。
分子ドッキングや MD シミュレーション前の前処理 として非常に役立ちます。
- Cl（塩素）、Na（ナトリウム）、K（カリウム）、Ca（カルシウム）、Mg（マグネシウム）の単独原子を削除
- 水分子（H₂O）を削除
- 処理後の分子を新しい SDF ファイルとして保存
- 現在のフォルダ内のすべての
.sdfファイルを自動処理
1. 除去する原子と水分子の定義
REMOVE_ATOMS = {"Cl", "Na", "K", "Ca", "Mg"}
REMOVE_SMARTS = ["[OH2]"] # H2O の SMARTS パターン
REMOVE_ATOMS：削除対象の 単独原子 を指定REMOVE_SMARTS：SMARTS（化学構造検索用の記法）で 水分子 (H₂O) を指定
2. Cl, Na, K, Ca, Mg の単独原子を削除
def remove_unwanted_atoms(mol):
"""Cl, Na, K, Ca, Mg の単独原子を削除"""
if mol is None:
return None
edit_mol = Chem.RWMol(mol) # 変更可能な分子オブジェクト
remove_idxs = [atom.GetIdx() for atom in mol.GetAtoms() if atom.GetSymbol() in REMOVE_ATOMS]
for idx in sorted(remove_idxs, reverse=True):
edit_mol.RemoveAtom(idx) # 指定した原子を削除
return edit_mol.GetMol() # RWMol → ROMol に変換
Chem.RWMol(mol)を使い、編集可能な分子オブジェクト に変換。- 分子内の全原子を走査し、
Cl,Na,K,Ca,Mgを持つ原子をリストアップ。 RemoveAtom()を使い、削除対象の原子を削除。
水分子 (H₂O) の削除
def remove_water_molecules(mol):
"""H2O の SMARTS パターンに一致する分子を削除"""
if mol is None:
return None
for smarts in REMOVE_SMARTS:
pattern = Chem.MolFromSmarts(smarts)
if mol.HasSubstructMatch(pattern):
return None # 水分子なら削除
return mol
- SMARTS で定義された 水分子 (H₂O) に一致する構造を持つ分子 を削除。
SDF ファイルの処理
def process_sdf(input_sdf, output_sdf):
"""SDF ファイルを処理し、不要な原子や水分子を削除"""
supplier = SDMolSupplier(input_sdf, removeHs=False) # SDF を読み込む
writer = SDWriter(output_sdf) # 出力用 SDF ファイルを作成
for mol in supplier:
if mol is not None:
mol = remove_unwanted_atoms(mol) # Cl, Na, K など削除
mol = remove_water_molecules(mol) # 水分子削除
if mol is not None:
try:
# **Kekulé 化のエラーを防ぐ**
Chem.SetAromaticity(mol, Chem.AromaticityModel.AROMATICITY_MDL)
Chem.SanitizeMol(mol, Chem.SanitizeFlags.SANITIZE_KEKULIZE)
# **ROMol に変換後、Kekulé 化をスキップして保存**
writer.write(mol)
except Chem.KekulizeException:
print(f"Warning: Failed to kekulize {input_sdf}. Skipping this molecule.")
writer.close()
print(f"Processed: {output_sdf}")
SDMolSupplierを使い、SDF ファイルを分子ごとにロード。remove_unwanted_atoms()で不要な原子を削除。remove_water_molecules()で水分子を削除。Chem.SanitizeMol(mol, Chem.SanitizeFlags.SANITIZE_KEKULIZE)を使い Kekulé 化のエラーを回避。SDWriter.write(mol)で 処理済み分子を新しい SDF ファイルに保存。
フォルダ内の .sdf ファイルを一括処理
# カレントディレクトリ内の全 .sdf ファイルを処理
input_sdf_files = glob.glob("*.sdf") # フォルダ内の .sdf ファイルを取得
output_sdf_files = [f"clean_{fname}" for fname in input_sdf_files]
for input_sdf, output_sdf in zip(input_sdf_files, output_sdf_files):
process_sdf(input_sdf, output_sdf)
glob.glob("*.sdf")で カレントディレクトリ内の全.sdfファイルをリスト化。clean_プレフィックスをつけた新しい.sdfファイルとして保存。- すべての
.sdfファイルを 一括処理 する。
水素、電荷付加、3D構造化、エネルギー最小化
以下で、低分子化合物の準備を順に行っていきます。状況に応じてclean_のついたsdfファイルの場所を変更してください。
###Add hydrogens
for i in *.sdf; do obabel $i -O $i -h; done
### Assign partial charges (Gasteiger charges)
for i in *.sdf; do obabel "$i" -O "$i" --partialcharge gasteiger; done
###Generate 3d
for i in *.sdf; do obabel $i -O $i --gen3d; done
#上記が通常だが、場合によっては時間がかかりすぎるものがあるので、以下に変更
ls *.sdf | parallel -j $(nproc) '
file="{}"
echo "Processing $file ..."
# ① 通常の3D変換（最大10分まで）
timeout 600 obabel "$file" -O "$file" --gen3d
if [ $? -eq 0 ]; then
echo "Standard --gen3d completed for $file"
exit 0
fi
echo "Timeout! Switching to --gen3d=fast for $file"
# ② 速い3D変換（最大10分まで）
timeout 600 obabel "$file" -O "$file" --gen3d=fast
if [ $? -eq 0 ]; then
echo "--gen3d=fast completed for $file"
exit 0
fi
echo "Timeout! Switching to 2D-to-3D approach for $file"
# ③ 2D→3D変換（最大10分まで）
obabel "$file" -O "2d_$file" --gen2d
timeout 600 obabel "2d_$file" -O "$file" --gen3d=fast
if [ $? -eq 0 ]; then
echo "2D to 3D conversion completed for $file"
exit 0
fi
echo "Timeout! Skipping stereochemistry processing for $file"
# ④ 立体化学情報をスキップ
obabel "$file" -O "$file" --gen3d=fast --ignore-stereo
echo "Final fallback used for $file"
'
###Then, Energy Minimize:
for i in *.sdf; obminimize -ff MMFF94 -n 1000 $i
1. 水素の付加
for i in *.sdf; do obabel "$i" -O "$i" -h; done
for i in *.sdf; do ...; doneは、カレントディレクトリ内の全てのSDFファイルに対してループ処理を行う構文。obabel "$i" -O "$i" -hで 水素原子を付加 し、元のSDFファイルを上書きする。hは 水素を追加するオプション で、適切なプロトネーション状態を考慮して水素原子を付与する。
この処理によって、通常 分子の電荷や結合数に基づいて水素が適切に追加 される。例えば、極性の高い酸素や窒素に適切な数の水素が付加される。
電荷の割り当て（Gasteiger電荷）
for i in *.sdf; do obabel "$i" -O "$i" --partialcharge gasteiger; done
-partialcharge gasteigerを使って Gasteiger電荷 を計算し、各原子に割り当てる。- Gasteiger電荷は、分子の電子密度分布を考慮し、経験的な手法で電荷を決定する方法の1つ。
この処理により、適切な電子分布が考慮された分子データが得られるため、後のエネルギー最適化がより物理的に意味のあるものになる。通常pdbqtファイルに変換する場合は自動的に電荷は割り当てられるが、今回はpdbqtファイルに変換しないので、この操作が必要。
3D座標の生成
for i in *.sdf; do obabel "$i" -O "$i" --gen3d; done
-gen3dは、2D構造しか持たない分子に3D座標を生成する。- Open Babel内部の ETKDG（Experimental Torsion Knowledge Distance Geometry）アルゴリズム を使用し、できるだけ実験値に近い立体構造を生成する。
この処理によって、SDFファイル内の 平面構造（2D SMILES 形式など）を、より現実的な3D構造へ変換 できる。
ls *.sdf | parallel -j $(nproc) '
file="{}"
echo "Processing $file ..."
# ① 通常の3D変換（最大10分まで）
timeout 600 obabel "$file" -O "$file" --gen3d
if [ $? -eq 0 ]; then
echo "Standard --gen3d completed for $file"
exit 0
fi
echo "Timeout! Switching to --gen3d=fast for $file"
# ② 速い3D変換（最大10分まで）
timeout 600 obabel "$file" -O "$file" --gen3d=fast
if [ $? -eq 0 ]; then
echo "--gen3d=fast completed for $file"
exit 0
fi
echo "Timeout! Switching to 2D-to-3D approach for $file"
# ③ 2D→3D変換（最大10分まで）
obabel "$file" -O "2d_$file" --gen2d
timeout 600 obabel "2d_$file" -O "$file" --gen3d=fast
if [ $? -eq 0 ]; then
echo "2D to 3D conversion completed for $file"
exit 0
fi
echo "Timeout! Skipping stereochemistry processing for $file"
# ④ 立体化学情報をスキップ
obabel "$file" -O "$file" --gen3d=fast --ignore-stereo
echo "Final fallback used for $file"
'
このスクリプトは Open Babel (obabel) を使用して .sdf 形式の分子データファイルを 3D 構造に変換するバッチ処理を行います。処理の流れとしては、以下の 4段階 で 3D 変換を試みます。
1. .sdf ファイルのリストを取得
ls *.sdf | parallel -j $(nproc) '
ls *.sdf：カレントディレクトリ内の.sdfファイルをリストアップparallel -j $(nproc) '...'：parallel：GNU Parallelを使用して並列処理を行うj $(nproc)：利用可能なCPUコア数(nproc)を指定して、並列で処理を実行
2. 通常の3D変換（最大10分）
timeout 600 obabel "$file" -O "$file" --gen3d
if [ $? -eq 0 ]; then
echo "Standard --gen3d completed for $file"
exit 0
fi
obabel "$file" -O "$file" --gen3d：-gen3dオプションを使用して 3D 構造を生成$file（入力ファイル）を上書き保存
timeout 600：- *最大600秒（10分）**で処理が終了しなければ自動終了
if [ $? -eq 0 ]; then exit 0; fi：obabelが正常終了 ($? -eq 0) すれば、そのまま終了- 失敗した場合は次のステップへ進む
3. 速い3D変換 (-gen3d=fast)
timeout 600 obabel "$file" -O "$file" --gen3d=fast
if [ $? -eq 0 ]; then
echo "--gen3d=fast completed for $file"
exit 0
fi
-gen3d=fastを使用して、高速な3D変換を試す- これも 10分のタイムアウト 付き
- 成功すれば終了、失敗すれば次の手段へ
4. 2D 構造を経由して 3D 変換
obabel "$file" -O "2d_$file" --gen2d
timeout 600 obabel "2d_$file" -O "$file" --gen3d=fast
if [ $? -eq 0 ]; then
echo "2D to 3D conversion completed for $file"
exit 0
fi
-gen2dを使って まず2D構造を生成- その後、2D構造 (
2d_$file) を 3D変換 - ここでも10分以内に変換が終わらなければ、さらに次の手段へ
5. 立体化学情報をスキップ
obabel "$file" -O "$file" --gen3d=fast --ignore-stereo
echo "Final fallback used for $file"
- 3D変換の最後の手段として、
-ignore-stereoを付けて 立体化学情報を無視して 変換 - これにより、 立体異性体の情報が欠落する可能性があるが、変換の成功率を高める
エネルギー最適化（エネルギー最小化）
for i in *.sdf; do obminimize -ff MMFF94 -n 1000 "$i"; done
obminimizeは、Open Babelの分子力場計算を用いた エネルギー最適化ツール。ff MMFF94は、MMFF94力場 を使用する指定。- MMFF94（Merck Molecular Force Field）は、有機分子のエネルギー計算に適した力場。
n 1000は、最大1000回の反復計算 を行い、エネルギー最小化を試みる。
この処理により、3D座標が物理的により妥当なコンフォメーションへ調整される。
以上で化合物ライブラリの準備は終わりです。
Grid boxの設定
Grid boxは実際にドッキングする場所を決める箱のことを指します。
MacではMGl toolsが利用できないので、UCFS chimeraを使って、Grid boxを決めます。
まずTools→General Controls→Command lineでCommand lineを有効にしてください。
次にSelect→Residueからリガンドの名前を確認してください。（ここではKIQ)
以下のような感じでコマンドラインにコードを書き、実行してください。
measure center :KIQ
リガンドの真ん中の座標がわかります。
続いて、Tools→Surface/Binding Analysis→Autodock vinaを押すと、次の画面が出てきます。
Resize search volume optionsにチェックを入れ、Centerに先ほど出力した数字、SizeにはGrid boxの大きさを指定し、ちょうどリガンドとその周辺が指定できるくらいの大きさに調節します。
これでGrid boxの設定ができました。CenterとSizeの値はメモってください。
Sminaによる分子ドッキング
ここからsminaをダウンロードします。smina.osx.12 というものがダウンロードできるので、こちらをApplicationフォルダに入れてください。以下のようにして、sminaコマンドの設定をしてください。
#Pathを通す
echo 'alias smina="/Applications/smina.osx.12"' >> ~/.zshrc
source ~/.zshrc
# 実行権限を付与
sudo chmod +x /Applications/smina.osx.12
# 確認
smina -h
このスクリプトは、smina という分子ドッキングツールのパスを設定し、実行権限を付与した後、動作確認を行うものです。
echo 'alias smina="/Applications/smina.osx.12"' >> ~/.zshrc
このコマンドは、~/.zshrc（zshの設定ファイル）に smina を簡単に実行できるようにエイリアスを追加 するもの。
sminaというコマンドを"/Applications/smina.osx.12"に置き換えて実行できるようになる。- 例えば、
smina -hと入力すると、実際には/Applications/smina.osx.12 -hが実行される。
source ~/.zshrc
設定をすぐに反映させるためのコマンド。
.zshrcを再読み込みすることで、エイリアスsminaが有効になる。- これを実行しない場合、ターミナルを再起動するまで
sminaは認識されない。
sudo chmod +x /Applications/smina.osx.12
smina.osx.12 に 実行権限（+x）を付与 する。
chmod +xは、ファイルを実行可能なプログラムとして認識させる。sudoを使うことで、管理者権限で実行し、システムのアプリケーションフォルダ内の変更を許可する。
smina -h
smina のヘルプメッセージを表示して、正常に動作するか確認 する。
- インストールが正しく完了し、エイリアスが機能していれば
smina -hでヘルプ情報が表示される。 - もし
command not found: sminaというエラーが出た場合、エイリアス設定が正しく反映されていない可能性があるので、source ~/.zshrcを再度実行するか、ターミナルを再起動する。
configファイルの設定
続いて、分子ドッキングに必要なconfigファイルの設定を行う。
config.txt
center_x = -18.53
center_y = -5.41
center_z = 59.97
size_x = 15
size_y = 15
size_z = 15
exhaustiveness = 8
num_modes = 9
energy_range = 3
内容は以下の通り。
| パラメータ | 説明 |
|---|---|
center_x, center_y, center_z | ドッキングボックスの中心座標（Å単位）事前に確認したリガンドの場所を設定する。 |
size_x, size_y, size_z | ドッキングボックスのサイズ（Å単位） |
exhaustiveness | 探索の強度（大きいほど精度向上、デフォルト=8） |
num_modes | 生成するポーズの数（デフォルト=9） |
energy_range | 受け入れるエネルギー範囲（kcal/mol, デフォルト=3） |
Sminaによるin silico screening
mkdir -p docking_results
#sminaの実行
for i in *.sdf; do
echo "Docking $i..."
smina -r 2i78_no_lig.pdb -l "$i" --config config.txt -o "docking_results/${i%.sdf}_docked.sdf" --log "docking_results/${i%.sdf}_log.txt" --scoring vina
done
このスクリプトは、Smina を用いて 分子ドッキング（ligand docking）を自動化 するためのものです。
カレントディレクトリ内の すべての .sdf ファイルを対象に、指定したタンパク質構造（2i78_no_lig.pdb）へドッキング を実行し、結果を docking_results フォルダに保存します。
mkdir -p docking_results
docking_results というフォルダを作成する。
pオプションは フォルダが既に存在する場合でもエラーを出さずに続行 できるようにするためのもの。- ドッキングの結果を保存するためのフォルダを用意する。
for i in *.sdf; do ... done
カレントディレクトリ内の すべての .sdf ファイル に対して、ループ処理を実行。
"$i"で各SDFファイルを順番に処理する。echo "Docking $i..."で現在処理中のファイルを表示（進捗確認用）。
smina -r 2i78_no_lig.pdb -l "$i" --config config.txt -o "docking_results/${i%.sdf}_docked.sdf" --log "docking_results/${i%.sdf}_log.txt" --scoring vina
Smina を用いて ドッキング計算を実行 する。
| オプション | 説明 |
|---|---|
-r 2i78_no_lig.pdb | 受容体（タンパク質） の構造ファイル |
-l "$i" | リガンド（化合物） のSDFファイル |
--config config.txt | 設定ファイル（ドッキングパラメータを指定） |
-o "docking_results/${i%.sdf}_docked.sdf" | ドッキング結果を保存するSDFファイル |
--log "docking_results/${i%.sdf}_log.txt" | ドッキングのログを保存 |
--scoring vina | AutoDock Vina のスコアリング関数を使用 |
${i%.sdf}_docked.sdf
${i%.sdf} は ファイル名から .sdf の拡張子を除いた部分 を取得する。
例えば ligand.sdf なら → ligand_docked.sdf というファイル名で保存される。
このスクリプトを実行すると、docking_results フォルダに 各リガンドのドッキング結果（座標）とスコア（ログ） が保存される。
すべての .sdf ファイルを一括処理できるため、大量のリガンドを自動でドッキング解析 するのに便利。
結果
結果は各_log.txtに保存してあります。これらの最も結合の高い（結合エネルギーの低い）もの順にソートを次のコードで行います。
import glob
import os
def parse_smina_log(log_file):
"""
Sminaのログファイルから mode1 の affinity 値を取得する。
"""
with open(log_file, 'r') as f:
for line in f:
if line.strip().startswith("1 "): # mode1 の行を取得
parts = line.split()
if len(parts) > 1:
try:
affinity = float(parts[1]) # Affinity値（kcal/mol）
return affinity
except ValueError:
pass
return None # mode1 のデータが見つからなかった場合
# Dockingログファイルを取得
log_files = glob.glob("*_log.txt")
# 各ログから affinity を取得
results = []
for log_file in log_files:
affinity = parse_smina_log(log_file)
if affinity is not None:
compound_name = log_file.split("/")[-1].replace("_log.txt", "")
results.append((compound_name, affinity))
# Affinity（結合エネルギー）の降順にソート（低いほど強い結合）
results.sort(key=lambda x: x[1])
# 結果をテキストファイルに出力
output_file = "docking_ranking.txt"
with open(output_file, "w", encoding="utf-8") as f:
f.write("ドッキング結果ランキング（結合力が強い順）\\n")
f.write("----------------------------------------\\n")
for rank, (compound, affinity) in enumerate(results, 1):
f.write(f"{rank}位: 化合物 {compound}, 結合エネルギー: {affinity:.2f} kcal/mol\\n")
# 結果を表示
print(f"{len(results)}個の化合物のドッキング結果をランキングしました")
print(f"結果は {output_file} に保存されました")
# トップ10の結果を表示
print("\\n=== トップ10の化合物 ===")
for rank, (compound, affinity) in enumerate(results[:10], 1):
print(f"{rank}位: 化合物 {compound}, 結合エネルギー: {affinity:.2f} kcal/mol")
以下のように出力されます。
6289個の化合物のドッキング結果をランキングしました
結果は docking_ranking.txt に保存されました
=== トップ10の化合物 ===
1位: 化合物 clean_drug1223, 結合エネルギー: -13.00 kcal/mol
2位: 化合物 clean_drug3004, 結合エネルギー: -12.60 kcal/mol
3位: 化合物 clean_drug2984, 結合エネルギー: -11.90 kcal/mol
4位: 化合物 clean_drug237, 結合エネルギー: -11.90 kcal/mol
5位: 化合物 clean_drug3074, 結合エネルギー: -11.80 kcal/mol
6位: 化合物 clean_drug538, 結合エネルギー: -11.70 kcal/mol
7位: 化合物 clean_drug92, 結合エネルギー: -11.70 kcal/mol
8位: 化合物 clean_drug4966, 結合エネルギー: -11.30 kcal/mol
9位: 化合物 clean_drug1884, 結合エネルギー: -11.20 kcal/mol
10位: 化合物 clean_drug214, 結合エネルギー: -11.10 kcal/mol
出力したclean_drug1223.sdf と2i78_no_lig.pdb をpymolなどで出力すると、分子ドッキングしていることがわかります。
※注意 こちら一番化合物が良かったclean_drug1223.sdf Obacunoneという化合物になりますが、pymolで判断したところ、立体構造の生成がうまくいっていません。ですので、再度立体構造別のツールでドッキングをするか、立体構造がきちんとしているリガンドかつ結合エネルギーが高いものを選ぶと良いと思います。
最後に
通常AutoDock vinaを使う場合は、MGltoolsを使う必要があり、このツールはMacでは使えません。そのため、MacにLinux環境を入れる必要があり、DockerでLinux環境を整えるということが一般的なようです。
M1 MacでDocker使ってAutodocktoolsを用いてドッキングシミュレーションをしよう – Qiita
しかし、私はうまくDockerでLinuxを構築する際、XQuartzがうまく実行できず、頓挫し、途方に暮れていました。そんな状況をXで呟いた時、@snowborderjack さんがsminaを教えてくれました。この場を借りて、感謝申し上げたい。ありがとうございました！
sminaは論文でも使われているので、今後も積極的に使っていきたい。
