オフラインでドッキングする話（T11パートBの備忘録）
T11パートBはオンラインWebサービスを用いたターゲットタンパク質に対するドッキングでした。残念ながら用いられていたSwissDock上のドッキング、OPAL共に現在では利用ができなくなっていました。そこで、趣旨とは異なってしまいますがオフラインでのドッキングを行って見たいと思います。
目標はこんな感じです。
| 用いるデータ | 概要 | 詳細 |
|---|---|---|
| 参照構造 | EGFR-リガンド共結晶構造 | PDB ID:3W33、但しKLIFSに登録されているデータ |
| ドッキングする化合物 | リガンドと構造類似の化合物群 | リガンド（ADP）をテンプレートとしてPubChemから取得 |
| やること | T11パートB | この記事 |
|---|---|---|
| 結合サイトの位置 | DoGSiteScorer | 同左 |
| タンパク質・リガンドの準備 | AutoDockTools | ODDT |
| ドッキング | OPAL上のAutoDockVina | PCにインストールしたAutoDockVina |
では、素人による適当解説やっていきますよ！
DogSiteScorerで良い感じのポケットを探す
Webサイトで遊ぶ
まずはドッキングを行うポケットの位置を決めます。ここはT11パートBでも実行できていたのでやり方を踏襲します。・・・が、やっぱりコードだけをみてもよくわからないのでWebサイトではどのような結果になるか見に行きたいと思います。
DoGSiteはProteinsPlusというサイトで提供されているソフトです*1。ProteinsPlusはハンブルク大学Computational Molecular Design Group（AMD、Matthias Rarey教授）で開発されたツールを利用できるように公開されているWebサーバーで、ポケットを探索するDoGSite以外にもタンパク質の構造を取り扱うためのツールが多数組み込まれています*2。
トップページ格好良い・・・
PDB IDでの検索も可能なようですが、今回はT11に従いKLIFSから取得したデータを使うため、Upload Proten(PDB format)にPDBファイルを投げ込みGo!。すぐに結果がでてこんなページが出てきます。
素敵・・・。
右側のパネルから利用できるサービスを選べるようになっています。DoGSiteScorerに行きましょう。色々とオプションがありますが今回はデフォルトのまま計算して見ます。
1−２分で結果が出ます。今回は10個のポケットが見つかりました。
ポケットのサイズに関する情報（体積、表面積）に加えてスコア（Drug Score、Simple Score）も出ています。スコアが高ければ医薬品の標的とするポケットとして良いのだと思います（論文を読んでいません。すみません）。確かにトップのポケット（P_0）は共結晶構造中のリガンドの位置と一致しています。
ポケットの情報はテーブルをスクロールした下部にあるDownloadから取得できます（zipファイル）。ありがたく頂戴して中身を見て見ましょう。
ダウンロードした情報の確認
ダウンロードしたディレクトリには以下のファイルが含まれています。
- 「parameter.txt」・・・計算の設定
- 「ジョブ名_desc.txtファイル」・・・ポケットの情報（Webサイト上に表示されているテーブル）。
- 「PocXlsDescriptors.txt」・・・ポケット情報の省略語の説明
- 「pockets」・・・各ポケットの座標等の情報のファイルを含むディレクトリ（.ccp4.gz、ポケットごとに別のファイル（今回は10個））
- 「residues」・・・各ポケットの残基等の情報のファイルを含むディレクトリ（.pdb、ポケットごとに別のファイル（今回は10個））
Webサイト上で画像表示されていたもののベースになる情報がそっくり取得できるようです。
ポケット情報のテキストファイル（2）はタブ区切りのテーブルのようでした。Pandasの出番だすな！
import pandas as pd df = pd.read_table('./DogSiteResults/Job_desc.txt') # ファイルの名前は長いので変えています print(df.shape) df.head() # (10, 50)
| name | lig_cov | poc_cov | lig_name | volume | enclosure | surface | depth | surf/vol | lid/hull | ... | MET | PHE | PRO | SER | THR | TRP | TYR | VAL | simpleScore | drugScore | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | P_0 | 0.0 | 0.0 | NaN | 1268.67 | 0.13 | 1422.18 | 26.50 | 1.121001 | - | ... | 2 | 3 | 2 | 1 | 2 | 1 | 0 | 4 | 0.64 | 0.817267 |
| 1 | P_1 | 0.0 | 0.0 | NaN | 371.84 | 0.17 | 681.06 | 13.45 | 1.831594 | - | ... | 2 | 1 | 3 | 0 | 1 | 0 | 1 | 2 | 0.23 | 0.642846 |
| 2 | P_2 | 0.0 | 0.0 | NaN | 338.56 | 0.07 | 524.19 | 13.65 | 1.548293 | - | ... | 1 | 1 | 1 | 3 | 1 | 2 | 1 | 2 | 0.18 | 0.625701 |
| 3 | P_3 | 0.0 | 0.0 | NaN | 234.94 | 0.10 | 473.53 | 10.68 | 2.015536 | - | ... | 0 | 2 | 2 | 1 | 0 | 0 | 1 | 1 | 0.11 | 0.421320 |
| 4 | P_4 | 0.0 | 0.0 | NaN | 207.36 | 0.12 | 353.27 | 11.60 | 1.703655 | - | ... | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 2 | 0.04 | 0.444417 |
5 rows × 50 columns
ポケットにどんな残基がいくつ含まれているかも含め、たくさんの情報が入っているようです。
具体的なポケットをPyMolで確認してみます。
ccp4ファイルはObjectMapCCp4として認識され、立方体が読み込まれました。ポケットを単位胞のように認識しているのでしょうか？・・・不勉強なのでわかりません。すみません。
pdbファイルのHEADERにはポケットの中心座標の情報も記載されています。例えば今回のポケットP_０では・・・
- HEADER Geometric pocket center at 0.63 18.01 37.64 with max radius 22.90
でした。この情報があればドッキング計算でポケットの位置を指定できそうです。
APIで実行
同じ作業をAPIで行ってみます。といってもT11パートBのコードを切りはりするだけですが・・・。
PDBファイルをアップロードするの複雑なのでそのままPDB IDを使ってやってみます（PDB ID: 3W33）。ProteinPlusはREST APIのドキュメンテーションもあり、DogSiteScorerはこちらです。POSTでジョブを作成し、GETで結果を取得します(json形式)。
import requests r = requests.post("https://proteins.plus/api/dogsite_rest", json={ "dogsite": { "pdbCode": "3W33", "analysisDetail": "1", "bindingSitePredictionGranularity": "1", "ligand": "", "chain": "" } }, headers= {'Content-type': 'application/json', 'Accept': 'application/json'} ) r.raise_for_status() # エラーに備える print(r.status_code) # 200
HTTP ステータス：200なのでうまくいっているようです。ドキュメンテーションには各ステータスコードの場合に得られる結果も記載されています。
実際にtext属性を確認してみると同じ結果となっています。
print(r.text) # {"status_code":200,"location":"https://proteins.plus/api/dogsite_rest/ジョブのID","message":"Job already exists"}
ジョブがうまく投げられたので計算結果を取得します。locationがジョブのIDなのでこちらを使います。
job_location = r.json()['location'] result = requests.get(job_location) print(result.status_code) # 200
こちらもうまくいったようです。結果は以下となっているそうです。
このうちresiduesの一番最初(インデックス0)が最もスコアの良いポケットだそうです。Webで結果をダウンロードした時と同じですね。
residuesはPDBファイル形式となっています。今欲しい情報は結合ポケットの中心の座標と半径で、こちらはHEADERのGeometric pocket center ~という行に含まれていました。こちらを取得します。
# PDBの取り出し best_pocket_pdb = requests.get(result.json()['residues'][0]).text # PDBから目的の情報を含む行を見つける # スペースで区切ったインデックス5,6,7（[5:8]）が中心の座標（x,y,z） # 最後（インデックス[-1]）が半径(radius) for line in best_pocket_pdb.splitlines(): line = line.strip() if line.startswith('HEADER') and 'Geometric pocket center at' in line: fields = line.split() center = [float(x) for x in fields[5:8]] radius = float(fields[-1]) print("center: ", center) print("radius: ", radius) # center: [20.37, 29.52, 10.77] # radius: 19.95
先ほどWebインターフェースで実行した場合と座標、半径が異なっています。KLIFSから取得したファイルとPDB IDを指定して行った違いのためでしょうか？同じ構造でも情報の取得減に応じて座標が異なる可能性があるというのは要注意そうです。
以降ではT11でKLIFSから取得したファイルを使うので、先にWebで得たポケットの情報を使うことにしたいと思います。
ドッキングの実行
ドッキングソフトの取得
結合ポケットの情報がわかったので実際のドッキングに進みたいと思います。ドッキングはAutoDockVinaで行いますが、AutoDockVinaを使うためにはタンパク質とリガンドをそれぞれ専用の形式（PDBQT）で準備する必要があります*3。T11パートAではAutoDockToolsを修正したものが使われていましたが、今回は別のツールODDT（Open Drug Dicovery Toolkit）を利用したいと思います。
ODDTは様々なオープンソースのツールの入出力方法を統合してもっと便利に使えるようにしましょう、という素敵なツールです。ツール多すぎて無理・・・って感じの私にはとてもありがたいです。詳細は以下をご参照ください。
- Wójcikowski, M., Zielenkiewicz, P., and Siedlecki, P. Open Drug Discovery Toolkit (ODDT): a new open-source player in the drug discovery field Journal of Cheminformatics 2015(7)26
- 文献URL
- ドキュメント
- GitHub
インストールはpipかcondaでできます。
# pipでインストールする場合 pip install oddt # condaでインストールする場合 conda install -c oddt oddt
AutoDockVina自体はこちらのページのDownloadから取得できます。計算化学.COMさんのこちらのページ にインストールと使用方法が詳しく書かれいるのでご参照ください*4。
タンパク質の準備
ではまずタンパク質だけを含むファイル（mol2）をODDTで読み込んでみます。ODDTではタンパク質の変数名をproteinとするようにとドキュメントに書かれていたので合わせておきます。*5
import oddt from oddt import toolkit protein = next(oddt.toolkit.readfile('mol2', "data/protein.mol2")) protein.protein = True type(protein) # oddt.toolkits.ob.Molecule
タンパク質の情報を読み込みオブジェクトの生成ができたので、こちらをpdbqtファイルに変換します。ODDTのdockingモジュールからAutodockVina.write_vina_pdbqt(molオブジェクト, 出力先のディレクトリ)を実行します。タンパク質の場合は引数flexible=Falseとしておく必要があります。
from oddt import docking oddt.docking.AutodockVina.write_vina_pdbqt(protein, "data", flexible=False) # 'data/3w33_A.pdbqt'
タンパク質の情報（PDB ID:3w33のChain A）を反映したファイル名で出力されました。変数proteinの中身をこちらに切り替えておきます。
protein = next(oddt.toolkit.readfile('pdbqt', "data/3w33_A.pdbqt")) protein.protein = True print(type(protein)) # oddt.toolkits.ob.Molecule
実際にはODDTはAutoDockVinaのドッキングを行う際にpdbqtへの変換も同時に行うので予め変換しておく必要はありませんが、ここでは使い方の練習ということでご容赦ください。
リガンドの準備
リガンドはT11パートAで取得したPubChem化合物を使います。SMILESで保存しているので、RDKitで３次元構造におこします。AllChem.EmbedMolecule()で３次元にしますが、性能が良いというETKDGv2()を指定して構造を生成します。*6
# SMILESファイルの読み込み SMILES_FILE = "data/similar_smiles.txt" with open(SMILES_FILE) as f: smiles = [line.strip() for line in f] # 必要なものをインポート from rdkit import Chem from rdkit.Chem import AllChem # ３次元構造を生成してリストに格納 mols = [] for s in smiles: mh = Chem.AddHs(Chem.MolFromSmiles(s)) AllChem.EmbedMolecule(mh, AllChem.ETKDGv2()) mols.append(mh) # 描画で確認 from rdkit.Chem import Draw Draw.MolsToGridImage(mols, molsPerRow=5)
立体感でてますね！SDFに書き出しておきます。
writer = Chem.SDWriter("data/similar_ligands.sdf") for mol in mols: writer.write(mol) writer.close()
ODDTでリガンドを読み込みPDBQTファイルに変換します。
# SDFの読み込み ligands = list(oddt.toolkit.readfile('sdf', 'data/similar_ligands.sdf')) # PDBQTに一つずつ変換 ligand_id = 0 for ligand in ligands: file_name = "ligand_" + str(ligand_id) oddt.docking.AutodockVina.write_vina_pdbqt(ligand, "data", name_id=file_name) ligand_id += 1
dataディレクトリの中に10個のリガンドpdbqtファイルができました。次のドッキングでは一番最初のリガンドを使って試そうと思うので読み込んでおきます。
ligand = next(oddt.toolkit.readfile('pdbqt', 'data/ligand_0.pdbqt'))
ドッキングの実行
結合ポケットの情報、タンパク質、リガンド構造の準備が完了したのでドッキングを実行します。ドッキングを行う空間はDoGSiteScorerのWeb版で実行しPyMolでも確認した以下のポケットです。
- 中心座標・・・ [0.63, 18.01, 37.64]
- サイズ・・・ 22.90 → 23にします。
from oddt.docking import AutodockVina # パラメーターの設定 DockingParams = AutodockVina.autodock_vina(protein=protein, size = (23, 23, 23), center = [0.63, 18.01, 37.64], exhaustiveness = 8, num_modes = 9, energy_range = 3, n_cpu = 3, executable='../../../../autodock_vina_1_1_2_mac/bin/vina')
引数executableはインストールしたAutodockVinaのlocalのpathを指定します。ここでは相対パスで指定しています。パラメーターの設定が完了したのでリガンドを渡してドッキングを実行します。
DockingResult = DockingParams.dock(ligand) len(DockingResult) # 9
9個の結果が生成されました。それぞれODDTのMolオブジェクトで、スコアはdataメソッドに格納されています。to_dict()で辞書に変換することができるので、辞書を経由してDataFrameに変換してみます。
data_list = [mol.data.to_dict() for mol in DockingResult] data_df = pd.DataFrame(data_list) data_df = data_df[['vina_affinity', 'vina_rmsd_lb', 'vina_rmsd_ub','vina_rmsd_input', 'vina_rmsd_input_min']] data_df
| vina_affinity | vina_rmsd_lb | vina_rmsd_ub | vina_rmsd_input | vina_rmsd_input_min | |
|---|---|---|---|---|---|
| 0 | -6.8 | 0.000 | 0.000 | 43.380367 | 43.380367 |
| 1 | -6.8 | 2.927 | 5.787 | 43.486763 | 43.486763 |
| 2 | -6.6 | 11.173 | 12.829 | 37.1132 | 37.1132 |
| 3 | -6.6 | 2.950 | 4.507 | 43.876957 | 43.876957 |
| 4 | -6.5 | 3.989 | 6.204 | 43.624863 | 43.624863 |
| 5 | -6.5 | 3.524 | 5.691 | 44.41204 | 44.41204 |
| 6 | -6.4 | 4.251 | 6.471 | 45.132126 | 45.132126 |
| 7 | -6.4 | 7.300 | 8.536 | 43.149113 | 43.149113 |
| 8 | -6.3 | 8.683 | 9.959 | 42.101715 | 42.101715 |
ここではvinaの出力に関係する列のみを取り出しています。vina_affinityとなっている列が親和性で最も良いもので-6.8となっています。DataFrameをcsvファイルに書き出しておきます。またドッキング結果のリガンドの情報もSDFに書き出しておきます。
# csv作成 data_df.to_csv("data/AutoDockVina_DockingResult.csv") # sdf作成 ODDTsdfWriter = oddt.toolkit.Outputfile('sdf', 'data/ODDT_DockingResult_ligand_0.sdf') for mol in DockingResult: ODDTsdfWriter.write(mol) ODDTsdfWriter.close()
引数のsdfの部分を変えればmol2等別の形式で出力することも可能です。
ODDTの出力におけるエラー
ODDTは非常に扱いやすいのですが、問題点としてAutoDockVinaでドッキング後のリガンドを出力すると原子、結合の帰属がおかしなことになっているということがあります。まず出力前の構造をJupyterNotebookで表示させたのちに、PDBQTで出力したリガンドをPyMolに表示させてみます。
DockingResult[0]
このように、骨格全体の座標としては保たれているものの、原子のタイプ、結合関係が壊れてしまっています。残念ながら出力のファイルタイプを変えてもこのエラーは修正されませんでした。PyMolで眺めて遊びたかったのに・・・
コマンドラインでAutoDockVinaを使う
ODDTで惜しいところまで行きましたが、私の能力では望みの構造の取得にまで到達できませんでした。残念・・・と思っていたらAutoDockVinaの使用方法のページ にコマンドラインで実行する方法が書いてありました。
コマンドはすごい苦手、というかさっぱりわからないのですが、ここまできたらやるしかない。
用意するものは３つ
- タンパク質のPDBQTファイル（ODDTで作成済み。flexible = Falseで作成していないと動かない）
- リガンドのPDBQTファイル（ODDTで作成済み）
- パラメーターを書いたconfigurationファイル(config.txt)
です。
config.txtの中身は以下のような感じ・・・
receptor = 3w33_A.pdbqt center_x = 0.63 center_y = 18.01 center_z = 37.64 size_x = 23 size_y = 23 size_z = 23 num_modes = 9
receptorにタンパク質のpdbqtファイル名、center_x,y,zはグリッドの中心座標、size_x,y,zはグリッドのサイズです。num_modesは最大いくつの結合ポーズを出すかを指定します。
ファイルの準備ができたら
- ターミナルでディレクトリに移動
- AutoDockVinaのインストール箇所へのPATHを通す。（
export PATH=$PATH:~~~~~） - コマンド実行！
vina --config conf.txt --ligand ligand_0.pdbqt --out out.pdbqt --log log.txt
と、出て結合ポーズを含むout.pdbqtとスコアなどの情報を含むlog.txtが出来ました！！早い！
PyMolで結果を眺めてみます。折角なので元の共結晶のリガンドを含むファイルと重ねてみます。
エモい！よくわからないけど！
スコアが良いものから悪いものへと順番に流れています。スコアの良いものは画面の右、悪いものは左側という印象でした。詳しい解析はT11パートCで扱われるはず・・・。
ちなみにAutoDockVinaは複数のリガンドをドッキングするためのシェルクリプトも用意してくださっています。こちらのページ からダウンロード出来ます。
chmodで実行権限を設定して実行すればまとめて処理できます。すごい便利ですね。食わず嫌いせずにシェルスクリプト勉強しないと・・・。
まとめ
以上、今回はドッキングを行ってみました。オンラインサービスは自分で環境を準備しなくてもよく、パワーのある計算機を使えるので素敵ですが、いつ止まるかわからないという不安があります。お手元のPCでも試せると便利かもしれません。心置きなく適当なことできますしね！
ODDTを導入したものの結果をうまく出力できなかったり、と相変わらずグダグダになってしまいましたが、T11パートBのオンラインではできなかったことをボンヤリとオフラインで実行できたので良し、ということにしていただければ幸いです。
色々と間違っていることがあると思うのでご指摘いただければ幸いです。
*1: DogSite文献 J. Chem. Inf. Model. (2010), 50, 11, 2041-2052
Bioinformatics (2012), 28, 15, 2074–2075
*3:J. Comput. Chem. 2010(31)455-61
*4:ドッキングシミュレーションのやり方【AutoDock vina】
*5: タンパク質ファイルの読み込みが遅く、セルが [＊] のまま進まないときはJupyter Notebookの拡張機能Variable Inspectorが原因の可能性があります。こちらの参考記事 【Jupyter Notebook】Jupyter notebookの動作が重い およびIssue と同じエラーかはわかりませんが、とりあえずnb_extensionsから外すと次のセルが実行できました。
