PyMOLでDocking Simulationできるのか -UNIXコマンド編-
SBDDの docking simulationに必要な探索範囲をboxとして可視化し、座標として算出する方法をこちらの記事で紹介しました。
この記事では可視化した座標を実際に探索するところまでやってみたいと思います。
フリーのdockingプログラムと言えばAutodock vinaが有名ですが、ここではより実行が簡単なvinaのfolkであるsminaを使いたいと思います。
GPUをサポートしたgninaも基本的な使い方は同じです。
ただし、これらはpythonではないので、そのままPyMOLで動かすことは出来ません。
Macならterminalからコマンドラインで動かす必要があります。
それをPyMOL上からやってしまおうと言うのが今回のお話です。
UNIXコマンド on PyMOL
PyMOLはコマンドも含めて基本的にpythonで動いています。
しかし、使い慣れた方は知ってると思いますが、cdやls、pwdなどのUNIXコマンドがいくつか使えます。
これは多分PyMOLのコマンドとして設定されているのだと思います。
試しにmkdirを入力してみても、下記のエラーが出ると思います。
mkdir
# PyMOL>mkdir
# NameError: name 'mkdir' is not defined
mkdirは定義されてないみたいです。
ついでにhelpコマンドも使ってみましょう。
help pwd
# PyMOL>help pwd
# DESCRIPTION
# Print current working directory.
# USAGE
# pwd
# SEE ALSO
# cd, ls, system
ちなみにhelpコマンドは自作したコマンドのhelpは出力しない仕様みたいです。
じゃあPyMOLでUNIXコマンドが動かないかと言うとそんなことなくて、osやsubprocessのライブラリをimportすることで使えるようになります。
jupyterとかでpythonを勉強してると時々出てくるライブラリです。
UNIXコマンドを作ってみる
docking simulationのような計算を実行すると、基本的に計算結果がどこかに出力されます。
計算結果の出力用にディレクトリがあるとデスクトップやディレクトリが散らからないので、PyMOLのコンソールからディレクトリを作成するmkdirのコマンドを作ります。
ついでにファイルを消去するrmコマンドも作ります。catのような複雑な動作のコマンドは難しそうなのでパスします。
# mkdir like unix command
def mkdir(dir:str):
os.makedirs(dir)
cmd.extend('mkdir', mkdir)
# rm command
def rm(file:str):
os.remove(file)
cmd.extend('rm', rm)
上記のスクリプトを保存してPyMOLコンソールからrunするか、シンプルにコピペしてコンソールに入力してもインストールできます。
試しに動かすとこんな感じ。
mkdir test
ls test
# test
testのディレクトリが出来てればOKです。
他にはsubprocessというライブラリでもUNIXコマンドを実行することができます。
こちらは複数の引数を取るようなコマンドでも実行できるようになります。
sminaはオプションが多いので、subprocessを使ってスクリプトを書いてます。
sminaを動かす
こちらの記事で紹介した探索範囲をboxで可視化するコマンドのスクリプトと統合したものがこちらです。
from pymol import cmd
from pymol.cgo import *
import itertools
from random import randint
import os
import subprocess
# mkdir like unix command
def mkdir(dir:str):
os.makedirs(dir)
cmd.extend('mkdir', mkdir)
# This script was described based on drawgridbox@PyMOL wiki. Users must check original script and license especially for commercial use.
# https://pymolwiki.org/index.php/DrawGridBox
"""
DESCRIPTION
Given selection, draw a box around it.
USAGE:
drawbox [selection, [padding, [lw, [r, [g, b]]]]]
PARAMETERS:
selection, the selection to enboxen
defaults to (all)
padding, defaults to 0
lw, line width
defaults to 2.0
r, red color component, valid range is [0.0, 1.0]
defaults to 1.0
g, green color component, valid range is [0.0, 1.0]
defaults to 1.0
b, blue color component, valid range is [0.0, 1.0]
defaults to 1.0
RETURNS
string, the name of the CGO box
NOTES
* This function creates a randomly named CGO box. The user can
specify the width of the lines, the padding and also the color.
"""
# draw gridbox
def grid(selection="enabled and organic", padding=4.0, lw=1.5, r=0.5, g=0.5, b=0.8):
([minX, minY, minZ], [maxX, maxY, maxZ]) = cmd.get_extent(selection)
print("Box dimensions (%.2f, %.2f, %.2f)" % (maxX - minX, maxY - minY, maxZ - minZ))
center_of_mass = [(minX + maxX) / 2, (minY + maxY) / 2, (minZ + maxZ) / 2]
print("Center of mass (%.2f, %.2f, %.2f)" % tuple(center_of_mass))
minX = minX - float(padding)
minY = minY - float(padding)
minZ = minZ - float(padding)
maxX = maxX + float(padding)
maxY = maxY + float(padding)
maxZ = maxZ + float(padding)
box = [
LINEWIDTH, float(lw),
BEGIN, LINES,
COLOR, float(r), float(g), float(b),
VERTEX, minX, minY, minZ, VERTEX, maxX, minY, minZ,
VERTEX, minX, maxY, minZ, VERTEX, maxX, maxY, minZ,
VERTEX, minX, minY, maxZ, VERTEX, maxX, minY, maxZ,
VERTEX, minX, maxY, maxZ, VERTEX, maxX, maxY, maxZ,
VERTEX, minX, minY, minZ, VERTEX, minX, maxY, minZ,
VERTEX, maxX, minY, minZ, VERTEX, maxX, maxY, minZ,
VERTEX, minX, minY, maxZ, VERTEX, minX, maxY, maxZ,
VERTEX, maxX, minY, maxZ, VERTEX, maxX, maxY, maxZ,
VERTEX, minX, minY, minZ, VERTEX, minX, minY, maxZ,
VERTEX, maxX, minY, minZ, VERTEX, maxX, minY, maxZ,
VERTEX, minX, maxY, minZ, VERTEX, minX, maxY, maxZ,
VERTEX, maxX, maxY, minZ, VERTEX, maxX, maxY, maxZ,
END
]
boxName = "box_" + str(randint(0, 10000))
while boxName in cmd.get_names():
boxName = "box_" + str(randint(0, 10000))
cmd.load_cgo(box, boxName)
return boxName
cmd.extend("grid", grid)
def search_space(selection="enabled and organic", padding=4.0):
([minX, minY, minZ], [maxX, maxY, maxZ]) = cmd.get_extent(selection)
print("Box dimensions (%.2f, %.2f, %.2f)" % (maxX - minX, maxY - minY, maxZ - minZ))
center_of_mass = [(minX + maxX) / 2, (minY + maxY) / 2, (minZ + maxZ) / 2]
print("Center of mass (%.2f, %.2f, %.2f)" % tuple(center_of_mass))
minX = minX - float(padding)
minY = minY - float(padding)
minZ = minZ - float(padding)
maxX = maxX + float(padding)
maxY = maxY + float(padding)
maxZ = maxZ + float(padding)
size_x = maxX - minX
size_y = maxY - minY
size_z = maxZ - minZ
center_x = (minX + maxX) / 2
center_y = (minY + maxY) / 2
center_z = (minZ + maxZ) / 2
return size_x, size_y, size_z, center_x, center_y, center_z
print("\nSmina Options:")
print("--center_x %.3f," % center_x, "--center_y %.3f," % center_y, "--center_z %.3f," % center_z,
"--size_x %.3f," % size_x, "--size_y %.3f," % size_y, "--size_z %.3f," % size_z)
cmd.extend("search_space", search_space)
def smina(ligand='enabled and organic', receptor='enabled and polymer.protein', grid=None, padding=4.0, n_pose=3, seed=0, comment=None):
# temporary file of receptor
rec='./rec.pdb'
# save pdb file of protein from enabled object
cmd.save(rec, receptor)
# temporary file of ligand
lig='./lig.pdb'
# save pdb file of protein from enabled object
cmd.save(lig, ligand)
if grid==None:
grid=ligand
else:
grid = grid + ' and polymer'
# defining grid box
([minX, minY, minZ], [maxX, maxY, maxZ]) = cmd.get_extent(grid)
center_of_mass = [(minX + maxX) / 2, (minY + maxY) / 2, (minZ + maxZ) / 2]
minX = minX - float(padding)
minY = minY - float(padding)
minZ = minZ - float(padding)
maxX = maxX + float(padding)
maxY = maxY + float(padding)
maxZ = maxZ + float(padding)
size_x = maxX - minX
size_y = maxY - minY
size_z = maxZ - minZ
center_x = (minX + maxX) / 2
center_y = (minY + maxY) / 2
center_z = (minZ + maxZ) / 2
if comment==None:
comment=''
print("\nExcuting Smina")
docking = [
'smina',
'-r', rec,
'-l', lig,
'-o', f'./docked_{comment}.sdf',
'--center_x', str(center_x),
'--center_y', str(center_y),
'--center_z', str(center_z),
'--size_x', str(size_x),
'--size_y', str(size_y),
'--size_z', str(size_z),
'--seed', str(seed),
'--num_modes', str(n_pose)
]
process = subprocess.Popen(docking, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = process.communicate()
print(stdout.decode())
print(stderr.decode())
cmd.load(f'./docked_{comment}.sdf')
cmd.extend("smina", smina)
from rdkit.Chem import PandasTools
def smina_score(docked='docked_.sdf'):
print(PandasTools.LoadSDF(docked))
cmd.extend('smina_score', smina_score)
いつも通りPyMOLのコンソールからスクリプトを実行します。
run ./grid.py
前回と同様に6NZPで動作を確認します。
fetch 6NZP
split_chain
disable 6NZP_B
grid
smina
下記のようにdockingの結果が読み込まれていれば成功です。
デフォルトでは表示しているタンパクの同じポケットでredockするようにスクリプトを書きました。
しかし、残念ながらtop rankのposeは元の構造と一致していませんでした。私の結果だと2番目のposeがオリジナルとよく一致しました。
sminaをbashなどで実行する場合、本来は下記のように記載します。
smina -r rec.pdb -l lig.pdb -o docked.sdf --autobox_ligand lig.pdb --seed 0
今回はPyMOLで処理できるように下記のようにpythonの引数を割り当てているので、以下に少し解説します。
def smina(ligand='enabled and organic', receptor='enabled and polymer.protein', grid=None, padding=4.0, n_pose=3, seed=0, comment=None):
ligand
ligand='enabled and organic'で、ligandはデフォルトで表示している有機分子を割り当てています。
他にもPyMOLのobjectとして読み込んでいる分子であれば、ligand=[object]でdockingすることができます。
rectptor
receptor='enabled and polymer.protein'で、receptorはデフォルトで表示しているタンパクを割り当てています。
ligandと同様にPyMOLで読み込んだobjectであればreceptorとして使えます。
注意として、apo体でない限りは下記のようにand polymer.proteinとする必要があります。
これは、receptorのファイルに残ったligandや水分子の位置を探索できないからです。
polymer.nucleicにすると核酸をreceptorとして残すこともできます。
smina receptor=6NZP and polymer.protein
grid
sminaでは探索する座標を設定する方法が大きく分けて2つあります。(詳しくは公式Documentを参照)
一つはautoboxでobjectの位置を探索させる方法、もう一つは探索範囲の座標を数値入力する方法です。
コマンドのオプションで言うと下記にあたります。--autobox_add argはgridコマンドpaddingにあたります。
--center_x arg
--center_y arg
--center_z arg
--size_x arg
--size_y arg
--size_z arg
--autobox_ligand arg
--autobox_add arg
詳細はスクリプトを確認してほしいですが、sminaコマンドのgrid=Noneオプションは表示しているオブジェクトの有機分子周辺を探索範囲にしています。
エチレングリコールのような添加剤が構造に含まれている場合はそちらも探索範囲に含まれるので、gridコマンドであらかじめ確認した方がいいです。
grid=のオプションで何かを指定する場合はand polymerとして処理されます。つまり特定のタンパク、アミノ酸、核酸の全体を指定します。
apo体でポケットを特定できている場合は、周辺のアミノ酸をseleとして選択して指定できます。
n_pose
n_pose=3で出力されるポーズ数をデフォルトで3にしています。
たくさん出力するほど十分に探索されますが、時間もかかります。
seed
seed=0は乱数値の指定です。指定することで同じdockingの結果が得られます。
ランダムな値を設定することでランダムな探索を実施します。その結果がいつも同様であれば、(結合情報あれば)信頼性は高そうだなという印象です。
(乱数値の設定忘れてました。思い出したら更新します。)
comment
comment=Noneで出力されるファイル(docked_.sdf)の末尾にコメントをつけます。
コメントをつけないとdocked_.sdfとしてdockingを実行するたびにファイルが更新されてしまうので注意です。
(overwriteの警告とかifで処理できるかも)
練習 -準備-
copilotによると、Serine/Threonin-Protein Kinase 2(SIK2)は構造情報がまだ報告されていない標的です(標的としての魅力はよく知りません)。 AlphaFold2のデータベースから予測構造をダウンロードし、ATPをdockingしてみましょう。 下記をクリックすると予測構造のpdbを直接ダウンロードできます(怖い人はしないこと)。 サンプルコードをTerminalで実行すると、SIK2のディレクトリを作ってpdbをダウンロードすることができます。
mkdir SIK2
cd SIK2
wget https://alphafold.ebi.ac.uk/files/AF-Q9H0K1-F1-model_v4.pdb
pymol AF-Q9H0K1-F1-model_v4.pdb
予測構造は比較的予測スコアが高めなようです。
下記のコマンドでpLDDTスコアを可視化します。
def plddt():
cmd.color('0x0053D6', 'enabled and b < 100')
cmd.color('0x65CBF3', 'enabled and b < 90')
cmd.color('0xFFDB13', 'enabled and b < 70')
cmd.color('0xFF7D45', 'enabled and b < 50')
util.cnc()
cmd.color('white', 'enabled and sc. and elem C')
cmd.extend('plddt',plddt)
plddt
予測構造は全長なので、機能にはおそらく不要かつ予測スコアの低い部分が大量に含まれています。
一旦非表示にして見やすくします。
hide cartoon, enabled and b < 50
show surface, enabled and b > 50
set surface_color, white
set transparency, 0.5
zoom ps. IEGTLGKGNFAVVKLGRHRITKTEVA
ちょっとわかりづらいんですが、kinaseのヒンジやP-loopが見えるようになると思います。
ここまでが準備とします。
練習 -AlphaFold2の予測構造にdocking-
まず下記の記事で紹介したPubChemから化学構造をダウンロードするコマンドをインストールしておきます。 相互作用解析用のコマンドもインストールして、後の解析に使います。
pc ATP
smina grid=AF-Q9H0K1-F1-model_v4 and b>50
fetch 1ATP
align 1ATP, AF-Q9H0K1-F1-model_v4
disable ATP
ATPをダウンロードしてdockingまで実行しました。
注意点として、予測構造そのままにdockingを試みると非ドメイン領域まで探索されて計算が終わりません。
そのため、grid=AF-Q9H0K1-F1-model_v4 and b>50として予測スコアの高い部分に限定して探索しています。
ついでにcAMP-dependent protein kinaseとATPの複合体(黄色)と重ね合わせています。
copilotによると、1ATPは世界で初めて報告されたkinaseとATPの複合体らしいです(ソース未確認)。
残念ながらATPのポケットにはdockingされませんでした！
気を取り直して、ATP pocketにdockingしてみましょう。
まずは探索範囲を可視化して確認しておきます。
これならATP pocketをちゃんと指定できそうです。
disable 1ATP
select sele, enabled and resi 30+95+100+160
grid sele, padding=0
disable docked_ box_*
enable ATP
smina grid=sele, padding=0, comment=local
line
enable 1ATP
hb
なんかいい感じにdockingされました！
ただし、細部を見てみると成功とはちょっと言えません。
今回はAF2の構造情報にMg2+が含まれていないため、三リン酸部分がポケット奥側に入ってしまっています。
コンフォメーションがおかしくなったせいか、アデニンがヒンジ部分とは相互作用位置に来ていませんでした。
AlphaFold2の構造はそのままdocking simulationに用いても上手くいくとは限らないと報告されています。
計算だけでそれっぽい見た目を得ることはできますが、妥当性を判断するドメイン知識が必要だとわかる実例になったかと思います。
おしまい。
Discussion