Zenn
🗳️

PyMOLでDocking Simulationできるのか -gridboxを作る編-

2024/12/17に公開
Python
pymol
docking
sbdd
tech

PyMOLはPythonベースのviewerです。　　
一方で、Schrödinger社のMaestroやCCGのMOEは有償ですが色んな計算をすることができる統合版のソフトウェアです。　　
計算まで実行できる無償版の統合ソフトウェアで言うとChimeraXもメジャーかなと思います。　　
最近ふと思い出したんですが、ChimeraXはAutodock Vinaを使ってdocking simulationまでサポートしてます。じゃあPyMOLでもやってみよう！というのが今回の内容です(なぜ。　　
ちなみに本記事ではAutodock Vinaのfolkで簡単に利用できるsminaを採用しています。　　
GPUをサポートしたgninaも書き換えてもらえば同じように動くと思います(未検証)。　　
それと、動作環境は基本的にM1 Mac miniです。　　
Linuxでも多分動きますが、Windowsのことは分かりません(無責任)。

introduction

低分子の分子シミュレーション、いわゆるSBDD(Structure-Based Drug Discovery)を行う場合(だけではありませんが)、事前の情報があるかどうかは成否に重要となります。　　
例えば、標的に対して結合するor活性のある分子が知られているか、構造情報が報告されているか、標的のどこに結合するか、などです。　　
創薬初期段階でこの情報が揃っていると、最適化までのスピードは格段に速くなります。　　
また、昨今では実際の物理的なライブラリの限界を突破する試みとして仮想分子ライブラリ(virtual library)のdockingを用いたin silicoスクリーニングの報告も増えてきているように思います。　　
2024年のノーベル化学賞の一つとなったAlphaFold2を代表とする構造予測技術と組み合わせて、大規模な探索をin silicoで実施し、wetでの検証を最小限に効率化させる取り組みは誰もが想像するところでしょう。　　
理想的には、in silicoスクリーニングと並行してアッセイ系を構築し、評価系が整った頃にはリード分子や候補分子が特定されているまで技術レベルを高めることだと思います。　　
話が若干逸れましたが、これらin silico創薬の実施に向けて前述の事前情報があると、wetでの検証項目が最小化され、確度とスピードが向上することが期待されます。

PyMOLでの課題の把握

冒頭述べた通り、PyMOLは基本的にviewerです。

  1. 標的タンパクの準備
  2. 探索範囲の指定 (本記事)
  3. ドッキングの実行

これらの計算はaddinで処理する必要があります。
残念ながら1の標的タンパクの準備で使いやすいプログラムを見つけられていないので、2と3について解決していきたいと思います。

gridboxを作る

docking simulationにおける探索範囲は大きく分けて、

  1. global search: 標的全体をscanning
  2. local search: 低分子が結合しそうなポケットに限定してscanning

の2つに分かれると思います。
どちらも一長一短で、前述の通りreferenceになる複合体構造情報があるならlocal searchすることが多いと思います。
標的のどこを探索するかを可視化し、指定するために、今回はPyMOL wikiに掲載されているdrawgridboxとcentroidを参考にしました。
https://pymolwiki.org/index.php/DrawGridBox
https://pymolwiki.org/index.php/Centroid
どちらもBSD(-2-Clause) Lisenceとのことなので注意してください。
前者はPyMOLの特定の座標を結んだgridboxを描画します。ちなみにgrid(格子)である必要がないので、スクリプトはシンプルなboxに書き換えてます。
後者は指定した座標の重心を計算します。
下記のスクリプトを保存したディレクトリでPyMOLを開き、スクリプトを実行します。

grid.py
from pymol import cmd
from pymol.cgo import *
import itertools
from random import randint

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)

# 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)
PyMOL
run ./grid.py

これでインストール完了です。早速、下記のスクリプトをPyMOLで実行して動かしてみましょう。

PyMOL
fetch 6NZP
split_chains
disable 6NZP_B
grid

6NZPは2023年に承認されたTYK2の阻害剤です。複合体構造情報が2019年に報告されています。
対称構造として報告されているので、split_chainsで分割した後にdisable 6NZP_BでB鎖を一旦消しています。　　
gridコマンドはデフォルトでは有機分子を認識しています。　　
また、この後にdockingをしたいのでpaddingで余白を少し取っています。padding=0のオプションでリガンドぴったりのboxを描けます。

PyMOL
grid padding=0

このboxの座標はPyMOLのコンソールに出力されています。search spaceコマンドでboxなしでも出力できます。

PyMOL
search_space

# PyMOL>search_space 
# Box dimensions (12.55, 10.08, 11.19)
# Center of mass (12.79, -3.08, 27.12)

既知リガンドがない標的の座標を指定するには

低分子の複合体が報告されている場合はgridコマンドのみで座標を特定できますが、実際にはapo体しか報告がないor得られていない標的もあると思います。　　
また、AlphaFold2のような予測構造では基本的には低分子の情報がありません。　　
仮に複合体構造情報がなくても結合の担保が取れた分子があればAlphaFold3で予測することも出来るかもしれませんが、それすらない場合はどうしましょう？　　
ある程度の規模のvirtual screeningでポケットにハマりそうな分子を選抜するなら、AlphaFold3よりも従来のdocking simulationの方が高速に処理できそうです。　　
(AF3をほとんど触ってないので実は詳細を知りません。有識者の皆さん、実際どうですか？)　　
仮に低分子が結合しそうなポケットがある場合は、周辺アミノ酸残基をseleとして選択してgridでboxを出力することができます。

PyMOL
grid sele

ポケットとかよくわからん、て言う場合はポケット予測のプログラムもいくつかあるので、そちらで絞り込むこともできると思います(今回は触れません)。　　
今回はdocking simulation向けに探索範囲を可視化・特定する方法を紹介したので、次回は実際にdocking simulationしてみたいと思います。

Discussion