こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル１１をもととしておりCC BY 4.0のライセンスに従っています。
トークトリアル１１は論文発表時のTeachOpenCADDには含まれておらず新しく追加されたものです。
４つのノートブックに分かれており、ここではその4つ目を訳出します。
パートCではパートBで行ったドッキング結果を解析しますが、パートBで使われていたオンラインサービスはすでに利用できなくなっていたため、ローカル環境ドッキングを行いました（実施方法はこちらの記事に記載いたしました）。
AutoDockVinaの２つの出力ファイルlog.txt、out.pdbqtをそれぞれファイル名を変更してlocal_vina.out、local_result.pdbqtとしてdataフォルダに格納し解析を行っています。
以下、訳。
トークトリアル 11 （パートC）
オンラインAPI/サービスを使った構造に基づくCADD
Developed at AG Volkamer, Charité
Dr. Jaime Rodríguez-Guerra, Dominique Sydow
このトークトリアルの目的
これは「オンラインWebサービス」についてのトークトリアルのパートCです。:
- 11a. キナーゼ阻害剤の候補をKLIFとPubChemで検索
- 11b. 11aで取得した候補化合物をターゲットタンパク質に対してドッキング
- 11c. 結果を評価し既知のデータと比較
入力構造を取得し、ドッキングが終わったので、結果が役に立つものかどうか評価します。
学習の目標
理論
実践
レファレンス
- タンパク質ーリガンド相互作用に関するレビュー(Int. J. Mol. Sci. (2016), 17, 144)
- PDBデータベースにおける非共有性結合相互作用の系統的な解析 (M. Med. Chem. Commun. (2017), 8, 1970-1981)
- タンパク質ーリガンド相互作用が医薬品の作用にとって鍵となることについての次の書籍の中の章（Klebe G. (eds) Drug Design. Springer, Berlin, Heidelberg. DOI: 10.1007/978-3-642-17907-5_4）
- Notebookのためのインタラクティブな可視化ソフト、NGLView (Bioinformatics (2018), 34, 1241–124)
- 自動的にタンパク質ーリガンド相互作用を見つけるためのPLIP(Nucl. Acids Res. (2015), 43, W1, W443-W447)
理論
タンパク質ーリガンド相互作用
リガンドの結合は、リガンドとタンパク質ポケットの表面との間、あるいはタンパク質ータンパク質境界面における非共有結合性相互作用によって主に支配されています。このプロセスは、静電的相補性と形状相補性、誘導適合（inducued fitting）、脱溶媒和過程、その他に応じて変化するものです。
文献からいくつか引用します。
José L. Medina-Franco, Oscar Méndez-Lucio, Karina Martinez-Mayorgaから改変して引用します。
タンパク質ーリガンド相互作用（protein-ligand interactions, PLIs）とタンパク質ータンパク質相互作用（protein-protein interactions、PPIs）を理解することは、分子の認識における核となる部分で、多くの科学分野において重要な役割を果たします。PLIsとPPIsは医薬品探索において幅広い分野で実践的な応用があります。含まれるものをいくつかあげれば、ドッキング、構造に基づくデザイン、フラグメント化合物と低分子化合物、そして他の種類の化合物のバーチャルスクリーニング、複合体のクラスタリング、そしてアクティビティクリフの構造的な解釈といったものがあり、他にもあります。
もちろん、これらの相互作用はいくつかの方法によって合理的に説明することが可能で、ドッキング結果の系統的な解析への道を開きます。例えばMed. Chem. Commun., 2017,8, 1970-1981改変して引用します。
私たちはPDBから、タンパク質と複合体を形成している低分子化合物で、分解能2.5Å以下の全てのX線構造を抽出し、その結果11,016個の複合体のコレクションを得ました。リガンドとして考慮するためには、化合物は低分子である、医薬品化学の応用にとって興味深いものである、といったいくつかの判断基準を満たしていなければなりません（緩衝液や結晶作成用の混合物の一部は除外します）。このコレクションには750,873のリガンドータンパク質の原子のペアが含まれており、ここでは原子のペアは4Å以内の距離にある２つの原子として定義されています。最も頻度の高いトップ100のリガンドータンパク質原子ペアは7つの相互作用のタイプに分類することができました（下の図を参照してください）。最も頻度高くみられた相互作用の中には、疎水性相互作用（hydophobic contact）、水素結合そしてπ-スタッキングといった、リガンドデザインにおいてよく知られており、広く使われている相互作用がありました。それらに続いて、弱い水素結合、塩橋、アミドスタッキング、そしてカチオンーπ相互作用といったものがありました。
タンパク質ーリガンド相互作用を自動的に評価するプログラムはいくつかあります。最も人気のあるプログラムの一つはPLIPで、Webサーバーが公に利用可能で無料で使えるPythonライブラリがあるというのが人気の理由です。文献に付随するsupporting informationはタンパク質ーリガンド相互作用を簡潔に解説しており、非常に簡単に理解できるようになっています。導入としては、次のパラグラフで十分でしょう。
PLIPはタンパク質残基とリガンドの非共有結合性相互作用をみつけるためのルールベースのシステムを使っています。タンパク質とリガンドの接触している原子間の非共有結合性相互作用の特徴を見つけるために、特定の相互作用に寄与することができる化学官能基（例、水素結合のドナーに必要とされるもの）と相互作用の幾何学的特徴（例、距離と角度の閾値）についての文献に基づく情報を使います。各結合サイトについて、このアルゴリズムは最初に、特定の相互作用のパートナーとなりうるタンパク質とリガンドの原子、あるいは原子団を探します。ついで、相互作用を形成するタンパク質とリガンドの条件に適合したグループに対して幾何学的なルールが適用されます。
PLIPのSIからのFigure。ここでは芳香族スタッキングが描かれています。
より詳細な情報が必要なら、PDFドキュメントをチェックしてください！
import requests
from IPython.display import HTML
pdf = "https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4489249/bin/supp_gkv315_nar-00254-web-b-2015-File003.pdf"
display(HTML(f"""
<iframe src="https://docs.google.com/viewer?url={pdf}&embedded=true"
frameborder="0" webkitallowfullscreen mozallowfullscreen
allowfullscreen width="900" height="600">
</iframe>
"""))
訳注(2020/05)
ここにはPLIPの文献PDFが表示されていますが、サイズの問題で反映させていません。PubMedのこちらのページでSUPPLEMENTARY DATAをご確認ください。
訳注ここまで
ドッキングソフトウェアは最終的な答えをあたえるものではありません！全ての結果は厳密に検証しなければなりません。例えば、このレビューL.Yang, J.Zhang, X.Che, Y.Q.Gao, in Methods in Enzymology, 2016 では次の指摘がなされています。
タンパク質ーリガンド相互作用の研究においては、リガンドの自由度も重要な役割を果たす可能性があります。しかしながら、リガンド/タンパク質相互作用の研究においては、通常非常に限られた数のコンフォメーションしか考慮されません。一方で、リガンド/基質ータンパク質の相互作用を調べ、予測する上で、リガンドのコンフォメーションの集合を徹底的に調べることが極めて重要であることを示す例が増えています。さらに、構造ベースの医薬品デザインにおいて、テストされたリガンドのコンフォメーションが続いて行う計算の基礎となります。特に、リガンドは作用するタンパク質に結合する前と結合した後で、非常に異なるコンフォメーションをとることがあり、多くのリガンドが局所最小エネルギーのコンフォメーションではタンパク質に結合していないことを示す証拠があります。そのような場合、自由なリガンド単体のコンフォメーションの不適切な解析は、誤った結合エネルギー予測につながり、さらに酵素反応が含まれている場合には問題となる可能性があります。したがって、自由なリガンド単体の溶液中の構造を徹底的にサンプリングすることが必要で、それは適切な強化されたサンプリング手法と組み合わせたMDシミュレーションによって費用対効果良く実現することができます。
これはありうる失敗する点のうちの一つです。
実践
結果の可視化
nglviewを使いましょう！Webベースの分子ビューワーでJupyter Notebook上で実行することができます！また、ボックスから取り出したPDBQTファイルとも互換性があります（が、最初のモデルだけしか読み込みません。。。これをどうするかは後ほどとりあげます）。
nglviewをインストールするには次を実行してください。
!pip install nglview ipywidgets>=7.5
!nglview install
!nglview enable
ノートブックでインタラクティブなGUIを作成するためにipywidgetsを使います。これにより、異なるリガンドをクリックし、クリックに応じてビューワーが更新されるようになります。特に、我々の小さなGUIには以下のことをさせたいと思います。
- Vinaの出力で報告されているように、ポーズと親和性のリストを表示する
- タンパク質構造をリボン表示、リガンドを Ball&Stick、そして周りの残基をlicorice(Stickのみ)で表示する
- 3Dの可視化がユーザーによるリストの異なるポーズの選択に応答するようにする
つまり、私たちは次のことをする必要があります。
- NGLビューワーが適切な描画をするよう呼び出す
- 結果のインタラクティブなテーブルを作成する（ヒント：
ipywidgets.Selectを使ってください）
- NGLビューワーとコミュニケーションできるイベントハンドラーを書く。つまりユーザーが新しいエントリをクリックすると、ディスプレイのリガンドと周囲の残基を更新するようにします（リボン表示は更新する必要がない）。
import pandas as pd
import time
import nglview as nv
from ipywidgets import AppLayout, Layout, Select, Button
Vinaの作成したPDBQTファイルはいくつかモデルを含んでいますが、nglviewは最初のものしか解析しません。この問題を回避する方法は単純です。ENDMDL行が見つかる度にスプリットすることでファイルを各モデルに分割します。
def split_pdbqt(path):
"""
複数のモデルのPDBQTを別々のファイルに分けます
"""
files = []
with open(path) as f:
lines = []
i = 0
for line in f:
lines.append(line)
if line.strip() == 'ENDMDL':
fn = f'data/results.{i}.pdbqt'
with open(fn, 'w') as o:
o.write(''.join(lines))
files.append(fn)
i += 1
lines = []
return files
Vinaの出力は結果のテーブルを含む単純なテキストファイルです。テーブルの解析は比較的簡単です。必要であれば、単純な可視化のためにpandas.DataFrameを返します。
def parse_output(out):
"""
Vina出力ファイルからDataFrameを作成します
"""
with open(out) as f:
data = []
for line in f:
if line.startswith('-----+'):
line = next(f)
while line.split()[0].isdigit():
index, *floats = line.split()
data.append([int(index)] + list(map(float, floats)))
line = next(f)
return pd.DataFrame.from_records(data,
columns=['Mode', 'Affinity (kcal/mol)', 'RMSD (l.b.)', 'RMSD (u.b.)'],
exclude=['Mode'])
Vinaの出力を解析
作成した関数をつかって次のテキストファイルを…
!cat data/local_vina.out
#################################################################
# If you used AutoDock Vina in your work, please cite: #
# #
# O. Trott, A. J. Olson, #
# AutoDock Vina: improving the speed and accuracy of docking #
# with a new scoring function, efficient optimization and #
# multithreading, Journal of Computational Chemistry 31 (2010) #
# 455-461 #
# #
# DOI 10.1002/jcc.21334 #
# #
# Please see http://vina.scripps.edu for more information. #
#################################################################
Detected 4 CPUs
Reading input ... done.
Setting up the scoring function ... done.
Analyzing the binding site ... done.
Using random seed: 267553816
Performing search ... done.
Refining results ... done.
mode | affinity | dist from best mode
| (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
1 -6.8 0.000 0.000
2 -6.8 5.072 7.461
3 -6.7 5.243 7.469
4 -6.6 4.035 6.317
5 -6.6 6.833 8.653
6 -6.6 8.876 10.836
7 -6.5 7.611 9.316
8 -6.5 1.094 2.109
9 -6.5 7.680 9.301
Writing output ... done.
綺麗にフォーマットされたpandas.DataFrameに変換することができます。
parse_output("data/local_vina.out")
|
Affinity (kcal/mol) |
RMSD (l.b.) |
RMSD (u.b.) |
| 0 |
-6.8 |
0.000 |
0.000 |
| 1 |
-6.8 |
5.072 |
7.461 |
| 2 |
-6.7 |
5.243 |
7.469 |
| 3 |
-6.6 |
4.035 |
6.317 |
| 4 |
-6.6 |
6.833 |
8.653 |
| 5 |
-6.6 |
8.876 |
10.836 |
| 6 |
-6.5 |
7.611 |
9.316 |
| 7 |
-6.5 |
1.094 |
2.109 |
| 8 |
-6.5 |
7.680 |
9.301 |
Vinaの出力には３つの列があります。
- 予測された結合親和性
- 対称性修正アルゴリズム（symmetry-corrected algorithm）を使った最良の解（親和性最小値）に対するRMSD
- 対称性修正（symmetry correction）を用いていない最良の解（親和性最小値）に対するRMSD
私たちは親和性にだけ興味があるので、次のセルではpandas.DataFrameのために列を一つだけ取得しています。
それでは、NGLビューワーのインスタンスを作成します。各タンパク質ーポーズのペアについて新しいインスタンスを作成する代わりに、同じ描画面を通して再利用し、必要なリガンドを隠すか表示します。最初に全てを読み込み、リガンドをそれぞれの親和性でラベルづけします。ビューワーに読み込まれた各分子を「component」とよびます。最初にタンパク質を読み込むので、component_0となります。リガンドはcomponent_1から始まり続いていきます。
import time
def create_viewer(protein, ligands, affinities):
"""
タンパク質と親和性のラベルをつけた全てのリガンドをもつnglview widgetを作成します
"""
viewer = nv.NGLWidget()
viewer.add_component(protein, ext="mol2")
label_kwargs = dict(labelType="text", sele="@0", showBackground=True, backgroundColor="black")
for ligand, affinity in zip(ligands, affinities):
ngl_ligand = viewer.add_component(ligand, ext="pdbqt")
ngl_ligand.add_label(labelText=[str(affinity)], **label_kwargs)
return viewer
そして最後に、次の下のセルで実際のGUIを作成します！
ipywidgets.AppLayoutレイアウトを使って水平に配置された２つのwidgetsからなります。
- セレクター（
ipywidgets.Select）
- NGLビューワーそのもの
ユーザーがセレクターで新しいエントリーをクリックすると、_on_selection_changeが呼ばれ以下のことを行います。
- 新しい値が前のものと異なるかを確認します。異なる場合、次に…
- 全てのリガンドを隠しビューをリセットします
- 選択したものを表示し、500msの格好いいアニメーションとともにカメラの中心に合わせます
- 新しいポーズの重心の5Å内にある残基のリストを更新するためにNGLビューワーにJavaScriptを実行します。
_RESIDUES_AROUND = """
var protein = this.stage.compList[0];
var ligand_center = this.stage.compList[{index}].structure.atomCenter();
var around = protein.structure.getAtomSetWithinPoint(ligand_center, {radius});
var around_complete = protein.structure.getAtomSetWithinGroup(around);
var last_repr = protein.reprList[protein.reprList.length-1];
protein.removeRepresentation(last_repr);
protein.addRepresentation("licorice", {{sele: around_complete.toSeleString()}});
"""
def show_docking(protein, ligands, vina_output):
ligands_files = split_pdbqt(ligands)
affinities = parse_output(vina_output)['Affinity (kcal/mol)']
viewer = create_viewer(protein, ligands_files, affinities)
selector = Select(options=[(f"#{i} {aff} kcal/mol", i) for (i, aff) in enumerate(affinities, 1)],
description="", rows=len(ligands_files), layout=Layout(width="auto"))
display(AppLayout(left_sidebar=selector, center=viewer, pane_widths=[1, 6, 1]))
def _on_selection_change(change):
if change['name'] == 'value' and (change['new'] != change['old']):
viewer.hide(list(range(1,len(ligands_files) + 1)))
component = getattr(viewer, f"component_{change['new']}")
component.show()
component.center(500)
viewer._execute_js_code(_RESIDUES_AROUND.format(index=change['new'], radius=5))
selector.observe(_on_selection_change)
_on_selection_change({'name': 'value', 'new': 1, 'old': None})
return viewer
上の関数ではPythonの応用をいくつか行なっています。もし興味があれば、下のヘッダーをクリックしてより詳細について見ることができます。
Python応用の説明
最初に、JavaScriptとPythonのコードを混ぜていることに気づいたかもしれません！これはNGLViewer.Viewer._execute_js_codeで提供されているように、ipywidgetsの基盤のおかげで可能になっています。ここで（JavaScriptのコードを含む）文字列を渡しておくと、後でNGL widgetのスコープで実行されます（thisはインタラクティブな描画面を差します）。パラメーター化するために、_on_selection_changeを呼ぶ度にフォーマット化されるテンプレート・プレースホルダーを追加しました。
この関数_on_selection_changeがユーザーのインタラクション（selectionボックス上でのクリック）をPythonの世界と繋げる糊の役割を果たします。selectionボックス上でクリックする度ごとに呼ばれ、引数changeを一つもちます。イベントのタイプと新旧両方の値を保持する辞書です（したがって比較し、それを取りあつかうことができます）。しかしながら、vewerを関数に含める必要もあります…そして、追加の値を渡すことはできません！
関数でviewerを使えるようにする一つの方法はshow_docking関数にその定義をネストすることです。そうすることで、Notebookの(global)scopeの中を探す必要なく外側の領域にアクセすることができるようになります。 このポストでscopeとclosureについて情報をもっと見ることができます。
ドッキング結果を見る
viewer = show_docking("data/protein.mol2", "data/local_result.pdbqt", "data/local_vina.out")
訳注(2020/05)
以下は画像ですが、ノートブック上ではインタラクティブなビューワーが表示されます。
訳注ここまで
viewer.render_image(),
(Image(value=b'', width='99%'),)
viewer._display_image()
import time
time.sleep(10)
タンパク質ーリガンド相互作用の解析
PLIPは自動化された解析のためのWebサーバー を提供していますが、残念ながらAPIはありません。標準のWeb UIを使っているかのようにHTML形式を使ってみることもできますが、ライブラリ自体がPython-3で使えるように準備されており、pipで簡単にインストールすることができるので、簡単のためローカルで使うことができます。
!conda install -y -c openbabel openbabel lxml
!pip install https://github.com/ssalentin/plip/archive/stable.zip --no-deps
from plip.modules.preparation import PDBComplex
from plip.modules.report import BindingSiteReport
from glob import glob
import os
import pandas as pd
PILIPはPDBファイルだけを受け入れますが、PDBQTファイルはPDBの拡張された記法なので、相互の変換は列とレコードのうち正しい部分を選択するだけの問題です。この場合、ATOM（タンパク質原子）とHETATM（リガンド原子）レコードにだけ興味があるので、追加のフィールド（66列以降）を廃棄します。ついで、単に連結を行うだけでそれらを合わせてまとめることができます。
def pdbqt_to_pdbblock(pdbqt):
lines = []
with open(pdbqt) as f:
for line in f:
if line[:6] in ('ATOM ', 'HETATM'):
lines.append(line[:67].strip())
return '\n'.join(lines)
def create_protein_ligand_pdbs(protein_pdbqt, ligands_pdbqt):
ligands_pdbqt = sorted(list(glob(ligands_pdbqt)), key=lambda s: int(s.split('.')[-2]))
protein_pdbblock = pdbqt_to_pdbblock(protein_pdbqt)
for i, ligand_pdbqt in enumerate(ligands_pdbqt, 1):
ligand_pdbblock = pdbqt_to_pdbblock(ligand_pdbqt)
with open(f"data/docked_protein_ligand.{i}.pdb", 'w') as f:
f.write(protein_pdbblock)
f.write("COMPND UNNAMED\n")
f.write(ligand_pdbblock)
ジャジャーン!できました！
create_protein_ligand_pdbs("data/protein.mol2.pdbqt", "data/results.*.pdbqt")
これでPDBファイルをPLIPに渡し、その魔法が使えるようになりました。BindingSiteReportクラスはPDBComplexの検出された各結合サイトを処理し、興味のあるフィールドをもつオブジェクトを作成します。
hydrophobichbondwaterbridgesaltbridgepistackingpicationhalogenmetal
これらのフィールドは（列名を保持する）<field>_featuresと（実際のレコードを保持する）<field>_infoに分けられます。オブジェクトを順に処理してgetattr()で正しい属性（attribute）の名前を取得すれば、概要をよくみるためにpandas.DataFrameに渡す辞書を構成することができます。以下をチェックしてください。
def analyze_interactions(pdbfile):
protlig = PDBComplex()
protlig.load_pdb(pdbfile)
for ligand in protlig.ligands:
protlig.characterize_complex(ligand)
sites = {}
for key, site in sorted(protlig.interaction_sets.items()):
binding_site = BindingSiteReport(site)
keys = "hydrophobic", "hbond", "waterbridge", "saltbridge", "pistacking", "pication", "halogen", "metal"
interactions = {k: [getattr(binding_site, k+"_features")] + getattr(binding_site, k+"_info") for k in keys}
sites[key] = interactions
return sites
interactions_by_site = analyze_interactions("data/docked_protein_ligand.1.pdb")
interactions_by_site
PLIP解析結果の辞書をみるにはここをクリックしてください
{'A:A:1': {'hydrophobic': [('RESNR',
'RESTYPE',
'RESCHAIN',
'RESNR_LIG',
'RESTYPE_LIG',
'RESCHAIN_LIG',
'DIST',
'LIGCARBONIDX',
'PROTCARBONIDX',
'LIGCOO',
'PROTCOO')],
'hbond': [('RESNR',
'RESTYPE',
'RESCHAIN',
'RESNR_LIG',
'RESTYPE_LIG',
'RESCHAIN_LIG',
'SIDECHAIN',
'DIST_H-A',
'DIST_D-A',
'DON_ANGLE',
'PROTISDON',
'DONORIDX',
'DONORTYPE',
'ACCEPTORIDX',
'ACCEPTORTYPE',
'LIGCOO',
'PROTCOO'),
(45,
'LYS',
'A',
1,
'A',
'A',
True,
'2.64',
'3.42',
'133.60',
True,
431,
'N3+',
2953,
'O3',
(3.232, 17.45, 38.708),
(4.782, 15.554, 36.326)),
(150,
'ASP',
'A',
1,
'A',
'A',
False,
'3.07',
'3.43',
'102.37',
True,
1470,
'Nam',
2955,
'O3',
(1.743, 18.722, 37.254),
(0.523, 19.864, 34.258)),
(45,
'LYS',
'A',
1,
'A',
'A',
False,
'2.66',
'3.10',
'106.87',
True,
422,
'Nam',
2951,
'O3',
(5.392, 21.015, 40.948),
(7.994, 19.469, 41.639)),
(72,
'LEU',
'A',
1,
'A',
'A',
False,
'3.65',
'4.01',
'103.83',
True,
694,
'Nam',
2945,
'Npl',
(3.257, 24.961, 34.128),
(2.841, 27.919, 36.804)),
(61,
'MET',
'A',
1,
'A',
'A',
False,
'2.75',
'3.76',
'161.43',
False,
2945,
'Npl',
595,
'O2',
(3.257, 24.961, 34.128),
(2.892, 28.611, 33.325)),
(149,
'THR',
'A',
1,
'A',
'A',
True,
'2.02',
'2.89',
'147.03',
False,
2955,
'O3',
1467,
'O3',
(1.743, 18.722, 37.254),
(0.131, 21.066, 36.729)),
(150,
'ASP',
'A',
1,
'A',
'A',
True,
'2.70',
'3.64',
'162.21',
False,
2953,
'O3',
1478,
'O3',
(3.232, 17.45, 38.708),
(1.632, 16.073, 35.738))],
'waterbridge': [('RESNR',
'RESTYPE',
'RESCHAIN',
'RESNR_LIG',
'RESTYPE_LIG',
'RESCHAIN_LIG',
'DIST_A-W',
'DIST_D-W',
'DON_ANGLE',
'WATER_ANGLE',
'PROTISDON',
'DONOR_IDX',
'DONORTYPE',
'ACCEPTOR_IDX',
'ACCEPTORTYPE',
'WATER_IDX',
'LIGCOO',
'PROTCOO',
'WATERCOO')],
'saltbridge': [('RESNR',
'RESTYPE',
'RESCHAIN',
'RESNR_LIG',
'RESTYPE_LIG',
'RESCHAIN_LIG',
'DIST',
'PROTISPOS',
'LIG_GROUP',
'LIG_IDX_LIST',
'LIGCOO',
'PROTCOO')],
'pistacking': [('RESNR',
'RESTYPE',
'RESCHAIN',
'RESNR_LIG',
'RESTYPE_LIG',
'RESCHAIN_LIG',
'CENTDIST',
'ANGLE',
'OFFSET',
'TYPE',
'LIG_IDX_LIST',
'LIGCOO',
'PROTCOO')],
'pication': [('RESNR',
'RESTYPE',
'RESCHAIN',
'RESNR_LIG',
'RESTYPE_LIG',
'RESCHAIN_LIG',
'DIST',
'OFFSET',
'PROTCHARGED',
'LIG_GROUP',
'LIG_IDX_LIST',
'LIGCOO',
'PROTCOO')],
'halogen': [('RESNR',
'RESTYPE',
'RESCHAIN',
'RESNR_LIG',
'RESTYPE_LIG',
'RESCHAIN_LIG',
'SIDECHAIN',
'DIST',
'DON_ANGLE',
'ACC_ANGLE',
'DON_IDX',
'DONORTYPE',
'ACC_IDX',
'ACCEPTORTYPE',
'LIGCOO',
'PROTCOO')],
'metal': [('RESNR',
'RESTYPE',
'RESCHAIN',
'RESNR_LIG',
'RESTYPE_LIG',
'RESCHAIN_LIG',
'METAL_IDX',
'METAL_TYPE',
'TARGET_IDX',
'TARGET_TYPE',
'COORDINATION',
'DIST',
'LOCATION',
'RMS',
'GEOMETRY',
'COMPLEXNUM',
'METALCOO',
'TARGETCOO')]}}
この辞書は２つのレベルからなります。最初のレベルは見つかった結合サイトです。一つのリガンドだけをドッキングしているので、一つだけ見つかることが予想されます。各結合サイトについて、８個のリストを含むもう一つの辞書があり、リストは各特定の相互作用のためのものです。各リストは最初の列に列名を持ち、つづいて（該当するものがあれば）データが格納されています。これを踏まえて、次のようにPandas.DataFrameを構築することができます。
df = pd.DataFrame.from_records(interactions_by_site["A:A:1"]["hbond"][1:],
columns=interactions_by_site["A:A:1"]["hbond"][0])
df
|
RESNR |
RESTYPE |
RESCHAIN |
RESNR_LIG |
RESTYPE_LIG |
RESCHAIN_LIG |
SIDECHAIN |
DIST_H-A |
DIST_D-A |
DON_ANGLE |
PROTISDON |
DONORIDX |
DONORTYPE |
ACCEPTORIDX |
ACCEPTORTYPE |
LIGCOO |
PROTCOO |
| 0 |
45 |
LYS |
A |
1 |
A |
A |
True |
2.64 |
3.42 |
133.60 |
True |
431 |
N3+ |
2953 |
O3 |
(3.232, 17.45, 38.708) |
(4.782, 15.554, 36.326) |
| 1 |
150 |
ASP |
A |
1 |
A |
A |
False |
3.07 |
3.43 |
102.37 |
True |
1470 |
Nam |
2955 |
O3 |
(1.743, 18.722, 37.254) |
(0.523, 19.864, 34.258) |
| 2 |
45 |
LYS |
A |
1 |
A |
A |
False |
2.66 |
3.10 |
106.87 |
True |
422 |
Nam |
2951 |
O3 |
(5.392, 21.015, 40.948) |
(7.994, 19.469, 41.639) |
| 3 |
72 |
LEU |
A |
1 |
A |
A |
False |
3.65 |
4.01 |
103.83 |
True |
694 |
Nam |
2945 |
Npl |
(3.257, 24.961, 34.128) |
(2.841, 27.919, 36.804) |
| 4 |
61 |
MET |
A |
1 |
A |
A |
False |
2.75 |
3.76 |
161.43 |
False |
2945 |
Npl |
595 |
O2 |
(3.257, 24.961, 34.128) |
(2.892, 28.611, 33.325) |
| 5 |
149 |
THR |
A |
1 |
A |
A |
True |
2.02 |
2.89 |
147.03 |
False |
2955 |
O3 |
1467 |
O3 |
(1.743, 18.722, 37.254) |
(0.131, 21.066, 36.729) |
| 6 |
150 |
ASP |
A |
1 |
A |
A |
True |
2.70 |
3.64 |
162.21 |
False |
2953 |
O3 |
1478 |
O3 |
(3.232, 17.45, 38.708) |
(1.632, 16.073, 35.738) |
さて、これらの相互作用をNGLビューワーで表しましょう。相互作用ポイント（pandas.DataFrameのLIGCOOとPROTCOO）の間に円柱（cylinder）を描き、color_mapに表されているように色分けします（RGBタプルを使います）。
color_map = {
'hydrophobic': [0.90, 0.10, 0.29],
'hbond': [0.26, 0.83, 0.96],
'waterbridge': [1.00, 0.88, 0.10],
'saltbridge': [0.67, 1.00, 0.76],
'pistacking': [0.75, 0.94, 0.27],
'pication': [0.27, 0.60, 0.56],
'halogen': [0.94, 0.20, 0.90],
'metal': [0.90, 0.75, 1.00],
}
def show_docking_and_interactions(protein, ligands, vina_output):
ligands_files = split_pdbqt(ligands)
affinities = parse_output(vina_output)['Affinity (kcal/mol)']
selector = Select(options=[(f"#{i} {aff} kcal/mol", i-1) for (i, aff) in enumerate(affinities, 1)],
description="", rows=len(ligands_files), layout=Layout(width="auto"))
app = AppLayout(left_sidebar=selector, center=None, pane_widths=[1, 6, 1], height="600px")
def _on_selection_change(change):
if change['name'] == 'value' and (change['new'] != change['old']):
if app.center is not None:
app.center.close()
app.center = viewer = nv.NGLWidget(height="600px", default=True, gui=True)
prot_component = viewer.add_component(protein, ext="pdbqt")
time.sleep(1)
lig_component = viewer.add_component(ligands_files[change['new']], ext="pdbqt")
time.sleep(1)
lig_component.center(duration=500)
interactions_by_site = analyze_interactions(f"data/docked_protein_ligand.{change['new']+1}.pdb")
interacting_residues = []
for site, interactions in interactions_by_site.items():
for interaction_type, interaction_list in interactions.items():
color = color_map[interaction_type]
if len(interaction_list) == 1:
continue
df_interactions = pd.DataFrame.from_records(interaction_list[1:], columns=interaction_list[0])
for _, interaction in df_interactions.iterrows():
name = interaction_type
viewer.shape.add_cylinder(interaction['LIGCOO'], interaction['PROTCOO'], color, [0.1], name)
interacting_residues.append(interaction['RESNR'])
res_sele = " or ".join([f"({r} and not _H)" for r in interacting_residues])
res_sele_nc = " or ".join([f"({r} and ((_O) or (_N) or (_S)))" for r in interacting_residues])
prot_component.add_ball_and_stick(sele=res_sele, colorScheme="chainindex", aspectRatio=1.5)
prot_component.add_ball_and_stick(sele=res_sele_nc, colorScheme="element", aspectRatio=1.5)
selector.observe(_on_selection_change)
_on_selection_change({'name': 'value', 'new': 1, 'old': None})
return app
show_docking_and_interactions はほとんどshow_docking関数とと同じですが、各選択イベントについてPLIP相互作用も計算します。この場合、イベントハンドラーはもう少しだけ関与しています。相互作用を示す円筒を手動で消すことができないので、各選択イベントについてNGLビューワーを作り直す必要があります。ですが、今回はJavaScriptは含まれていません。純粋なPythonだけです。ですが、ビューワーのwidgetを初めから作り直し、分子を読み込み、相互作用を計算するにはより時間が必要となり、体感としてはそれほどダイナミックなものにはなりません。
ドッキングと相互作用の表示
app = show_docking_and_interactions("data/protein.mol2.pdbqt", "data/local_result.pdbqt", "data/local_vina.out")
app
訳注(2020/05)
以下は画像ですが、ノートブック上ではインタラクティブなビューワーが表示されます。
訳注ここまで
app.center.render_image(),
(Image(value=b'', width='99%'),)
app.center._display_image()
ディスカッション
このトークトリアルの最後のパートでは、Vinaによるドッキングの結果を可視化し解析する方法を学びました。他のパートと異なり、このNotebookはローカル環境で実行され、Webサーバーは呼びませんでした。しかしながら、とてもたくさんのWebテクノロジーを取り扱っており、主にNGLビューワーとipywidgetsレイアウトに基づいています。
はじめに、選択したタンパク質ーリガンドポーズとその周囲でインタラクティブにビューワーを更新するために(observe()で登録されているように）イベントハンドラーを使いました。次に、ビューワーで相互作用を描くためにPLIPの自動化された解析を使って可視化を拡張しました。
クイズ
- 観察したドッキング解析において最も重要だったのはどのような相互作用でしたか？文献との一致は見られましたか？またそれはなぜでしょうか？
- 疎水性相互作用と水素結合の間の主な違いは何でしょうか？またどの点で類似しているでしょうか？
- パートBで学んだ内容を応用してPLIPのWebサーバーを使う関数を描いてください。メインの
formエントリーとactionのエンドポイント、リクエストを投げるのに必要なinputタグの場所を見つけるためにHTMLソースコードをチェックする必要があります。
訳以上
追記
以上です。今回も色々と誤りがあると思います。ご指摘いただければ幸いです。
以下、今回のトークトリアルでよく分からなかった箇所とエラーについて調べた内容を記載します。
- 追記 1) AutoDockVinaの結果、RMSDのsymmetry correction に関して
RMSDについてsymmetry correctionを行ったか否かの２種類の記述があります。以下は、私の理解です。
リガンドの各原子を区別してラベルしている場合、分子に対称性があると、化学的には同じ構造でも原子のラベルとしては異なるという可能性があります（ex. ベンゼンの６つの炭素原子にラベルをつける場合）。ドッキングの結果の比較として、同じリガンドの異なる結合ポーズをRMSDで評価しようとする場合、各原子のラベルの対応を元に計算できます。ですが、対称性がある分子では、実際には重なりのよい（RMSDの小さい）ポーズでも、ラベルの付け方が複数通り可能なために原子同士が離れている（RMSDが大きくなる）と計算してしまいます。この問題を考慮するか否か、というのがsymmetry correctionを行うか否か、という違いだと思います。
PLIPでPDBを読み込む際に下記のエラーが出ました。
~/.pyenv/versions/anaconda3-5.3.1/envs/teachopencadd/lib/python3.6/site-packages/plip/modules/supplemental.py in read_pdb(pdbfname, as_string) in read_pdb(pdbfname, as_string)
396 if os.name != 'nt': # Resource module not available for Windows
397 maxsize = resource.getrlimit(resource.RLIMIT_STACK)[-1]
--> 398 resource.setrlimit(resource.RLIMIT_STACK, (min(2 ** 28, maxsize), maxsize))
399 sys.setrecursionlimit(10 ** 5) # increase Python recoursion limit
400 return readmol(pdbfname, as_string=as_string)
ValueError: current limit exceeds maximum limit
使用するシステムのリソースを設定する箇所でのエラーです。こちらはどうもPLIPの問題ではなく、Mac OSとPython 3.6の組み合わせで生じる問題のようです。Pythonにissue（ 課題34602:python3 resource.setrlimit strange behaviour under macOS ） が出ていました。PLIPの該当箇所をコメントアウトしGitHubの修正履歴を参考に以下のように書き換えました。
def read_pdb(pdbfname, as_string=False):
"""Reads a given PDB file and returns a Pybel Molecule."""
pybel.ob.obErrorLog.StopLogging() # Suppress all OpenBabel warnings
"""エラー箇所
if os.name != 'nt': # Resource module not available for Windows
maxsize = resource.getrlimit(resource.RLIMIT_STACK)[-1]
resource.setrlimit(resource.RLIMIT_STACK, (min(2 ** 28, maxsize), maxsize))
"""
# 変更ここから
if sys.platform == 'darwin':
try:
import resource
except ImportError:
pass
else:
soft, hard = resource.getrlimit(resource.RLIMIT_STACK)
newsoft = min(hard, max(soft, 1024*2048))
resource.setrlimit(resource.RLIMIT_STACK, (newsoft, hard))
# 変更ここまで
sys.setrecursionlimit(10 ** 5) # increase Python recoursion limit
return readmol(pdbfname, as_string=as_string)
これでエラーが出なくなりました。
ビューワーを画像にしてしまっているためわかりにくいですが、ドッキング結果のリガンドは結合におかしな点がみられます。ドッキング結果を出力したPDBファイルには原子の座標のみで、結合関係の情報が入っていないようでしたので、その辺りが影響しているのかもしれません。
終わりに
これでトークトリアル１１含め全てのトークトリアル終了となります。長々とお付き合いいただきありがとうございました！GitHubにノート入れておきましたのでよろしければご利用ください。
github.com
私は、CADDに関連する項目についてはぼんやりと用語を聞いたことがある程度だったので、各トピックごとに理論と実践の両方が扱われていて、まとめて勉強する良い機会となりました。一つの標的EGFRに絞って、色々なトピックでアプローチすることで、様々な視点から情報を捉えるというのは面白いやり方ですね。恐らく内容としては初歩の初歩で、CADDのほんの入り口に過ぎないのだと思いますが、全体の流れと道しるべ的な専門用語を拾えたので、何が話題になっているか程度には理解できるようになっているといいなと思います。ではでは！
