magattacaのブログ

日付以外誤報

TeachOpenCADD トピック10 〜結合サイトの類似性とオフターゲット予測〜

こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリトークトリアル10をもととしておりCC BY 4.0のライセンスに従っています。
GitHubはてなMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!

以下、訳。

トークトリアル 10

結合サイトの類似性とオフターゲット予測

Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin

Angelika Szengel, Marvis Sydow and Dominique Sydow

Note: このノートブックはセルを順番に実行してください。一度に全てのセルを実行することもできますが、PyMolのイメージが思った通りにならない可能性があります。

このトークとリアルの目的

このトークトリアルでは、タンパク質全体および結合サイトの構造的類似性を使ってオフターゲットを予測します。オフターゲットとは意図せず医薬品のターゲットとなっているタンパク質で、望まぬ副作用につながったり、あるいは、当初の意図とは異なる望ましい医薬品の適応に繋がることがあります(ドラッグリポジショニング)。結合サイトの比較の主要なステップについて議論し、基本的な手法、つまり構造間の幾何学的な差異(2つの構造の平均二乗偏差(root mean square deviation, rmsd))を実装します。

学習の目標

理論

  • オフターゲットタンパク質
  • コンピューターによるオフターゲット予測:結合サイト比較
  • 類似性の簡単な指標としてのペアワイズ平均二乗偏差(Pairwise RMSD)
  • チロシンキナーゼ阻害剤、イマチニブ、

実践

  • 着目するリガンドを読み込み、可視化する(イマチニブ/STI)
  • PDBから全てのタンパク質-STI複合体を取得
  • PDB構造を可視化
  • PDBの構造をアラインメント(タンパク質全体)
  • ペアワイズRMSDの取得(タンパク質全体)
  • PDBの構造をアラインメント(結合サイト)
  • ペアワイズRMSDの取得(結合サイト)

レファレンス

結合サイト比較:

イマチニブ(Imatinib):

理論

オフターゲットタンパク質

オフターゲットは、医薬品あるいは医薬品の代謝物(の一つ)と相互作用するタンパク質で、ターゲットタンパク質として意図されていないものです。オフターゲットによって起こされる分子の応答は、望まぬ副作用につながる可能性があり、あまり害の無いものから非常に有害な副作用まであります。オフターゲットが生じる主な理由は、本来のターゲット(on-target)とオフターゲット(off-target)が互いの結合サイトに類似の構造的モチーフを共に有しており、それゆえ類似のリガンドに結合することができる、ということです。

コンピューターによるオフターゲット予測:結合サイト比較

コンピューター支援による潜在的なオフターゲットの予測の目的は、潜在的な危険性をもつ医療用物質を開発するリスクを最小にすることです。結合サイトの類似性評価にはいくつかのアルゴリズム的アプローチがありますが、常に次の3つの主要なステップで構成されています。

  1. 結合サイトの符号化(エンコーディング:結合サイトを様々な記述子テクニックでエンコードし、ターゲットデータベースに格納する
  2. 結合サイトの比較:様々な類似性指標を使って、クエリの結合サイトをターゲットデータベースと比較する
  3. ターゲットランキング:適切なスコア化に基づきターゲットを順位づけする

様々な類似性指標と既存のツールについてのより詳細な情報は、結合サイト比較に関する次の素晴らしいレビュー2本を参照して下さい(Curr. Comput. Aided Drug Des. (2008), 4, 209-20 そしてJ. Med. Chem. (2016), 9, 4121-51)。

f:id:magattaca:20200506005325p:plain

Figure 1: 結合サイト比較手法の主な工程(Dominique Sydowによる図)。

類似性の簡単な指標としてのペアワイズ平均二乗偏差

類似性のスコア化のための簡単かつ直接的な手法は平均二乗偏差(RMSD)の計算値を使うことです。RMSDはアラインメントをとった2つの構造の原子間の距離の2乗の平均の平方根となります (Wikipedia: RMSD)。

2つの構造間で比較するそれぞれの原子を見つけるために、まず、配列に基づく方法あるいは配列とは独立なアラインメントアルゴリズムによって、2つの構造のアラインメントを取る必要があります (Book chapter: Algorithms, Applications, and Challenges of Protein Structure Alignment in Advances in Protein Chemistry and Structural Biology (2014), 94, 121-75)。

チロシンキナーゼ阻害剤、イマチニブ

キナーゼはATPからタンパク質へとリン酸を移す酵素で、それによりシグナル伝達や代謝、タンパク質の制御といった様々な細胞のプロセスをコントロールします。 これらのキナーゼが(ゲノムの変異のせいで)恒常的に活性化していると、制御プロセスを歪ませ癌の原因となることがあります。抗がん剤の例としてイマチニブがあげられます(Nat. Rev. Clin. Oncol. (2016), 13, 431-46)。癌の治療に使われる低分子チロシンキナーゼ阻害剤で、より具体的には慢性骨髄性白血病(Chronic Myeloid Leukaemia、CML))と消化管間質腫瘍(GastroIntestinal Stromal Tumour、GIST)の治療に使われます。

イマチニブの選択性は完全なものでは無く、その主要なターゲット以外のチロシンキナーゼも標的とすることが示されています。このことはドラッグリポジショニングに使われており、つまり、イマチニブは異なる癌種に対しても承認されています (J. Biol. (2009), 8, 10.1186/jbiol134)。しかしながら、アレルギー反応や感染、出血、頭痛といった望まぬ副作用も示します (Drugs.com: Imatinib)。

実践

以下では、イマチニブに結合するPDB構造を取得しフィルタリングします。イマチニブの結合したタンパク質の(構造の解かれたタンパク質との)構造類似性を調べます。類似性の指標は(簡単な類似性指標として)ペアワイズRMSD計算を使います。潜在的なオフターゲットをみつけるための最初の評価として、このシンプルな手法が利用できることを示します。

import os
import pprint
import pickle
import glob
import time

from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw, AllChem
from rdkit.Chem import PyMol
IPythonConsole.ipython_useSVG=True

import pandas as pd
import numpy as np
import random
import matplotlib.pyplot as plt
import seaborn as sns

from pypdb import *
from biopandas.pdb import PandasPdb

対象のリガンドを読み込み可視化する(イマチニブ(Imatinib/STI))

イマチニブ(よく使われる略語: STI)のSMILES形式はChEMBLデータベース(ChEMBL: Imatinib) あるいはPDBデータベース(PDB: STI)から、略語「STI」で取得することができます。ここでは単純に手作業でChemical Component Summaryテーブルの「Isomeric SMILES」エントリーから文字列をコピーし、リガンドを読み込みます。

訳注(2020/05)
イマチニブに関するPDBのChemical Component Summarのリンクはこちら(RCSB PDB: STI )です。
訳注ここまで

sti = Chem.MolFromSmiles('CN1CCN(Cc2ccc(cc2)C(=O)Nc2ccc(C)c(Nc3nccc(n3)-c3cccnc3)c2)CC1')
Draw.MolToImage(sti)

f:id:magattaca:20200506005407p:plain

STIの3次元構造を調べるために、オープンソースツールのPyMolを使います(トークトリアルT8のイントロダクションを参照してください)。PyMolでSTIを眺める前に、3次元座標を計算する必要があります。

まず、化合物に水素原子を追加しますが、水素原子はSMILES形式でいつも明示的に記述されているとは限りません。次に、距離幾何学(ディスタンスジオメトリー法、Distance Geometry)を使って化合物の初期座標を取得し、UFF(Universal Force Field)力場を使って化合物構造の最適化を行います。

sti_mol = Chem.AddHs(sti)
sti_mol

f:id:magattaca:20200506005426p:plain

AllChem.EmbedMolecule(sti_mol)
AllChem.UFFOptimizeMolecule(sti_mol)
sti_mol

f:id:magattaca:20200506005445p:plain

これで、PyMolに進む準備ができました!

ターミナルからPyMolを開きます(ここではJupyter notebookからコントロールしています)。

os.popen('pymol -R')

objPMV = PyMol.MolViewer()経由でPyMolとJupyter notebookを接続するために、PyMolが完全に立ち上がるまで待つ必要があります。

# エラーへの対処方法:PyMolが読み込まれるまでまつ

nrTry = 0  # 現在の試行回数
ttw = 10  # 待ち時間(秒単位)

while nrTry < ttw:  # PyMolが読み込まれオブジェクトが保存されるまで試行する
    nrTry += 1
    try:  
        objPMV = PyMol.MolViewer()  # PyMolオブジェクトの保存
        break  # PyMolが読み込まれたらループを止める        
    except ConnectionRefusedError:  # PyMolがまだ読み込まれていない場合の例外処理
        time.sleep(1)  # 待つ...
    
if nrTry == ttw:  # ttw試行の後:エラーメッセージを表示
    print("Error: PyMol did not start correctly.\n" +
          "Try again and/or check if PyMol is installed completely.")

STIをPyMolに読み込み、このJupyter notebookにPyMolフレームを表示させます。

# PyMolを初期化(すでにPyMolにオブジェクトを読み込んでいるといけないので初期化する)
objPMV.server.do("reinitialize")

# PyMolに化合物を読み込み表示
objPMV.ShowMol(sti_mol)
# sticksで表示
objPMV.server.do('cmd.show("sticks","molecule")')
    
# PyMolの背景を白色に設定
objPMV.server.do("bg_color white")
# 現在のPyMolフレームのレイトレースイメージを作成
objPMV.server.do("ray 900, 900")

objPMV.GetPNG(h=300)

f:id:magattaca:20200506005532p:plain

PDBから全てのタンパク質-STI複合体を取得

PDBのような公開データベースでイマチニブ/STIを調べ、ターゲットタンパク質として報告されているタンパク質を検索することができます。PDBでは、イマチニブは通常STIとして省略されています。以下では、両方の名称を検索し結果を統合します。

PDBを検索

まず、対象のリガンド(STI)に結合するタンパク質をPDBから全て取得します。

search_dict = make_query('STI')  # PDBでリガンドSTIに結合するタンパク質を検索
found_pbd_ids = do_search(search_dict)

print(found_pbd_ids)
print("\nNumber of structures connected with STI in the PDB: " + str(len(found_pbd_ids)))
#   ['1AVU', '1AVW', '1AVX', '1BA7', '1FPU', '1IEP', '1M52', '1OPJ', '1R8N', '1R8O', '1T45', '1T46', '1XBB', '2BEA', '2BEB', '2DRE', '2ESU', '2ET2', '2HYY', '2OIQ', '2PL0', '3EA7', '3EA8', '3FW1', '3GVU', '3HEC', '3I2A', '3I2X', '3K5V', '3M3V', '3MS9', '3MSS', '3OEZ', '3PYY', '3QQM', '3S8J', '3S8K', '4BKJ', '4CSV', '4GCN', '4GCO', '4H9W', '4HA2', '4R7I', '4TLP', '5FNW', '5FNX', '5FZU', '5FZY', '5FZZ', '5G00', '5MQT', '6HD4', '6HD6', '6I0I', '6JOL', '6KTN', '6NPE', '6NPU', '6NPV', '6OTU', '6WE5']
#  Number of structures connected with STI in the PDB: 62

検索対象のリガンドに使う名称(ここではイマチニブ(imatinib))によって検索結果が変わることがあることに注意してください。

search_dict2 = make_query('Imatinib')  # PDBでリガンド Imatinib結合するタンパク質を検索
found_pbd_ids2 = do_search(search_dict2)

print(found_pbd_ids2)
print("\nNumber of structures connected with Imatinib in the PDB: " + str(len(found_pbd_ids2)))

# ['1IEP', '1M52', '1OPJ', '1T46', '1XBB', '2F4J', '2GQG', '2HYY', '2OIQ', '2PL0', '2XYN', '3EL7', '3EL8', '3FW1', '3G0E', '3G0F', '3G6G', '3G6H', '3GVU', '3HEC', '3HEG', '3K5V', '3MS9', '3MSS', '3OEZ', '3OF0', '3PYY', '3QLF', '3QLG', '4BKJ', '4CSV', '4R7I', '5MQT', '6HD4', '6HD6', '6JOL', '6KTM', '6KTN', '6NPE', '6NPU', '6NPV']
    
#  Number of structures connected with Imatinib in the PDB: 41

両方の検索結果を統合し、重複しない(ユニークな)エントリーだけを残します。

pdb_ids = list(set(found_pbd_ids + found_pbd_ids))
print("Number of structures connected with STI/Imatinib in the PDB: " + str(len(pdb_ids)))
#  Number of structures connected with STI/Imatinib in the PDB: 62

訳注(2020/05)
set()で2つのリストの和集合をとっています(重複を削除)。set型となるため再度list()でリスト型にしています。
訳注ここまで

PDBデータセットをフィルタリング

以下のクライテリアに従ってデータセットをフィルタリングするために、pypdbget_entity_info関数を使ってPDB構造のメタ情報を取得します。

  1. 実験手法でフィルタリング (xray)
  2. 分解能でフィルタリング (3 Å以下)
  3. (簡単にするため)1本鎖のPDB構造のみを残す(single chain)
  4. Imatinibに結合した構造のみを残す(例、Imatinibに関連づけられているが、結合はしていないPDB構造が除かれる)
  5. 2019年以前に登録されたPDB IDのみを残す(トークトリアル公開時点のリソースにデータセットをあわせる)

PDBを検索する方法についての情報はトークトリアル 8を参照してください。

# PDBからメタ情報を取得
entity_info = []
for i in pdb_ids:
    entity_info.append(get_entity_info(i))
entity_info[0]
{'Method': {'@name': 'xray'},
 'Entity': {'@id': '1',
  '@type': 'protein',
  'Chain': [{'@id': 'A'}, {'@id': 'B'}]},
 'structureId': '6I0I',
 'bioAssemblies': '1',
 'release_date': 'Wed Sep 04 00:00:00 PDT 2019',
 'resolution': '1.45'}
# リストをデータフレームに変換
entity_info_pd = pd.DataFrame(entity_info)
a = [int(i.split(" ")[5]) for i in entity_info_pd["release_date"].tolist()]
entity_info_pd.head()
Method Entity structureId bioAssemblies release_date resolution
0 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': [{'@... 6I0I 1 Wed Sep 04 00:00:00 PDT 2019 1.45
1 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': {'@i... 5FZU 1 Wed Jul 06 00:00:00 PDT 2016 2.43
2 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': {'@i... 2PL0 1 Tue Oct 09 00:00:00 PDT 2007 2.80
3 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': {'@i... 1R8N 1 Tue May 25 00:00:00 PDT 2004 1.75
4 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': [{'@... 3M3V 1 Wed Mar 23 00:00:00 PDT 2011 2.70
# 1. 実験手法でフィルタリング
entity_info_pd = entity_info_pd[entity_info_pd["Method"] == {'@name': 'xray'}]

# 2. 分解能でフィルタリング
entity_info_pd = entity_info_pd[entity_info_pd["resolution"].astype(float) <= 3.0]

# 3. 1本鎖の構造のみを保持(簡単化のため)
entity_info_pd = entity_info_pd[[type(i) == dict for i in entity_info_pd["Entity"]]]

entity_info_pd = entity_info_pd[[type(i["Chain"]) == dict for i in entity_info_pd["Entity"]]]

pdb_ids = entity_info_pd["structureId"].tolist()

print("Number of structures after filtering: " + str(len(pdb_ids)))
# Number of structures after filtering: 24
entity_info_pd.head()
Method Entity structureId bioAssemblies release_date resolution
1 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': {'@i... 5FZU 1 Wed Jul 06 00:00:00 PDT 2016 2.43
2 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': {'@i... 2PL0 1 Tue Oct 09 00:00:00 PDT 2007 2.80
3 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': {'@i... 1R8N 1 Tue May 25 00:00:00 PDT 2004 1.75
5 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': {'@i... 1T45 1 Tue Jun 15 00:00:00 PDT 2004 1.90
6 {'@name': 'xray'} {'@id': '1', '@type': 'protein', 'Chain': {'@i... 1AVU 1 Wed Oct 28 00:00:00 PST 1998 2.30

訳注
Chainに関する情報はカラム「Entity」の辞書型の要素のうち、「Chain」をkeyとするvalueとして格納されています。
Chainの数に応じて、
* Chainが1つの時 → 「'Chain': {'@id': 'A'}」 → valueのtypeはdict * Chainが2つの時 → 「'Chain': [{'@id': 'A'}, {'@id': 'B'}]」 → valueのtypeはlist

のようになっており、複数のChainがある場合はvalueのtypeがlist(辞書を要素とするリスト)になります。
従って、フィルタリング3行程目では「type(i["Chain"]) == dict」 とすることで「Chainが1つのもののみ」を残していると思われます。
訳注終わり

以下では、BioPandasというパッケージを使用します。BioPandasにはpandasデータフレームに生物学的高分子(PDBからmol2ファイルまで)の分子構造を読み込むために役立つ関数があります。PDBファイルを使った作業を容易にするためにPandasPdbオブジェクトを使います。

# 4. Imatinibの結合した構造だけを保持

def check_if_ligand_present(pdb_id, ligand_name):
    ppdb = PandasPdb().fetch_pdb(pdb_id)  # PDB(原子情報、座標)を取得
    return sum(ppdb.df["HETATM"]["residue_name"] == ligand_name) > 0  # エントリー STIがあるかどうかb確認

entity_info_pd = entity_info_pd[[check_if_ligand_present(pdb_id, "STI") for pdb_id in pdb_ids]]  # 関数を適用

pdb_ids = entity_info_pd["structureId"].tolist()

print("Number of structures after filtering: " + str(len(pdb_ids)))
#   Number of structures after filtering: 9

訳注(2020/05)
BioPandasのドキュメンテーションこちらです。
PandasPdbオブジェクトのアトリビュートPnadasPdb.dfDataFrameの辞書となっており、PDBファイルにpandasデータフレームのようにアクセスできます(BioPandas Tutorial)。ここではPDBファイルのHETATMレコードにアクセスし残基名(residue_name)でリガンドの名称を取得しています。

また、check_if_ligand_presentは引数ligand_nameが一つでもあればTrue、なければFalseを返す関数です。上記のセルではリスト内包表記で「STI有り無しのboolリスト」を取得し、これを使ってTrueとなる行のみをDataFrameから抽出しています。
訳注ここまで

# 5. 2019年以前に登録されたPDB DIのみを保持

entity_info_pd = entity_info_pd[[int(i.split()[5]) < 2019 for i in entity_info_pd["release_date"].tolist()]]

pdb_ids = entity_info_pd["structureId"].tolist()

print("Number of structures after filtering: " + str(len(pdb_ids)))
#  Number of structures after filtering: 8

訳注(2020/05)
「release date」は「例:Wed Jul 06 00:00:00 PDT 2016」という形式なので、半角スペースでsplitした6番目(インデックス5)が「年」になります。
訳注ここまで

# 6. 手作業で目視で調べたあとで3GVUを除く(STIを2つ含む:以下の自動化ワークフローに適していない)
pdb_ids.remove("3GVU")
pdb_ids
# random.shuffle(pdb_ids)  # IDの順序を変えたい時に使用してください

# ['2PL0', '3FW1', '3HEC', '4R7I', '1XBB', '4CSV', '1T46']

フィルタリングしたPDB IDの保存

さらに解析を進めるため、フィルタリングしたデータセットPDB IDを保存します(後ほど、このファイルに従ってPDB IDを処理するPyMolのPythonスクリプトを使います)。

pickle.dump(pdb_ids, open("../data/T10/pdb_ids.p", "wb"))

PDB構造の可視化

まず、タンパク質データセットの3次元構造を目視で調べるためPyMolに全構造を読み込みます。

このJupyter notebook上でのPyMolフレームの固定化したイメージ形式での可視化に加えて、PyMolアプリケーション上で直接3次元構造を眺め、インタラクティブに検証することをお勧めします。PyMolアプリケーションはこのトークトリアルと並行して立ち上り、操作されているはずです。

# PyMolの再初期化
objPMV.server.do("reinitialize")

# タンパク質ファイルの読み込み
for pdb_id in pdb_ids:
    cmd = "fetch " + pdb_id
    objPMV.server.do(cmd)

# タンパク質をcartoonとして表示(使用しているpymolのversionに応じて必要な可能性があります)
objPMV.server.do('cmd.show("cartoon","all")')

# オブジェクトを隠す
objPMV.server.do("hide polar_contacts")

# 背景を白色に設定
objPMV.server.do("bg_color white")
# 水分子とイオンを削除
objPMV.server.do("remove solvent")
# 中心に置き、ズーム
objPMV.server.do("center all")
objPMV.server.do("zoom all")
objPMV.server.do("ray 400,400")

# Jupyter noteboomにPyMolフレームを表示
objPMV.GetPNG(h=500)

png

このイメージは美しくカラフルでくるくると巻いていますが、まだ有益な情報を与えてくれません。次のステップではこの構造を互いに重ね合わせます(アラインメント)。

PDBの構造をアラインメント(タンパク質全体)

PyMolはさまざまなレベルの配列類似性と構造類似性に適した、異なるアラインメント手法を提供しています:

  • superコマンドはかなりの構造類似性(decent structural similarity)のあるタンパク質に適しています(配列には依存しない)。
  • alignコマンドはかなりの配列類似性(decent sequence similarity)のあるタンパク質に適しています(配列に依存)。
  • cealign コマンドはほとんどあるいは全く配列類似性がない(little to no sequence similarity)タンパク質に対して頑健な手法です(配列依存性は中間)

このトークトリアルではalignコマンドを選択し、(配列に基づいて)RMSDを最小化することで構造を重ね合わせます。

注:このアプローチは類似の配列を持つ構造に有利なように解析されます(alignコマンドはかなりの配列類似性をもつタンパク質のペアでパフォーマンスがより良くなります)。配列類似性が低いタンパク質の比較にはsuperあるいはcealignコマンドの方が選択肢として良いかもしれません。(タンパク質ペアの配列や構造類似性がわからない)自動化されたワークフローにの場合、3つの全ての手法に基づいてRMSDを計算し、最も結果の良いものをさらなる解析のために残すことが解決策となるかもしれません。

まず、pdb_idsリストの全ての構造を一番最初の構造に対してアラインメントをとります。

immobile_pdb_id = pdb_ids[0]
for mobile_pdb_id in pdb_ids[1:]:
    objPMV.server.do("align " + mobile_pdb_id + ", " + immobile_pdb_id)

objPMV.server.do("center all")
objPMV.server.do("zoom all")
objPMV.server.do("ray 400,400")

objPMV.GetPNG(h=500)

f:id:magattaca:20200506005836p:plain

タンパク質のうちの一つ、3FW1は他のタンパク質と比較して重なりが悪くなっています。重なりの良いタンパク質を表示させるために、このタンパク質を隠します。

objPMV.HideObject("3FW1")

objPMV.server.do("center all")
objPMV.server.do("zoom all")
objPMV.server.do("ray 400,400")

objPMV.GetPNG(h=500)

f:id:magattaca:20200506005855p:plain

多くの部分(例、ヘリックス)では構造の重なりの程度は高いですが、他の部分では重なりは低いか悪くなっています。

ペアワイズRMSDの取得(タンパク質全体)

(PyMolとPyMol Pythonスクリプトでは可能ですが)このJupyter notebook上ではPyMolからRMSD値を取得する関数を見つけられなかったので、PyMolのPythonスクリプトを呼びます。PyMol Pythonスクリプトでは全てのタンパク質のペアについて重ね合わせと最終的なRMSDのリファインを行います。結果としてえらえれるRMSD値をさらなる解析のために保存しておきます。

まず、このタスクのために使う、PyMolのPythonスクリプトを見ます。

f = open("./pymol_scripts/pymol_align_proteins.py")
file_content = f.readlines()
for i in file_content:
    print(i, end="")
"""
    '''
    This script aligns all pairs of input proteins.
    '''
    # 訳者実行環境の都合上pymolのpathを追加するために以下2行追記
    import sys
    sys.path.append('/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages')
    # 追記終わり。以降オリジナルスクリプト
    
    import pymol
    import pickle
    import pandas as pd
    import sys
    import numpy as np
    
    # Launch PyMol
    pymol.finish_launching()
    
    # Load PDB IDs
    pdb_ids = pickle.load(open("../data/T10/pdb_ids.p", "rb"))
    
    # Create DataFrame to store align results
    align_df = pd.DataFrame(index=pdb_ids, columns=pdb_ids)
    
    def align_structures(m1, m2, align_df):
      # Fetch PDB structures
      pymol.cmd.fetch(m1, "i")
      pymol.cmd.fetch(m2, "m")
    
      # Align proteins
      aln = pymol.cmd.align(mobile="m", target="i", quiet=1, object="{}_{}".format(m1, m2))
    
      # Save align results to DataFrame
      align_df.loc[m1, m2] = aln
    
      # Save sequence alignment to file
      pymol.cmd.save(filename="../data/T10/alignment/alignment_proteins_" + "{}_{}".format(m1, m2) + ".aln", selection="{}_{}".format(m1, m2))
    
      # Reinitialize PyMol for next structure pair to be loaded
      pymol.cmd.reinitialize()
    
      return align_df
    
    # Iterate over all structure pairs
    for n, immobile in enumerate(pdb_ids):
      for mobile in pdb_ids[n:]:
          align_df = align_structures(mobile, immobile, align_df)
          align_df = align_structures(immobile, mobile, align_df)
    
    # Quit PyMol
    pymol.cmd.quit()
    
    # Save results to file
    pickle.dump(align_df, open("../data/T10/align_df_proteins.p", "wb"))
"""
os.popen("python ./pymol_scripts/pymol_align_proteins.py")

訳注
アラインメント結果が保存されていない場合、Pythonスクリプトのエラーの可能性があります。
私の実行環境ではターミナル上でModuleNotFoundError: No module named 'pymol' が確認できたため、PyMolのPathを追加したところ実行することができました。
訳注終わり

このPyMol Pythonスクリプトは全てのペアワイズRMSDの結果をファイルに保存します。それでは結果をJupyter notebookに読み込みます。

align_df_proteins = pickle.load(open("../data/T10/align_df_proteins.p", "rb"))

align_df_proteinsは各ペアワイズの比較についてPyMolのalignコマンドの戻り値、つまり7項目のタプル、を含むデータフレームです:

  1. リファイン後のRMSD
  2. リファイン後の重ね合わされた原子の数
  3. リファインメントサイクルの試行数
  4. リファイン前のRMSD
  5. リファイン前の重ね合わされた原子の数
  6. アラインメントスコアの生の値
  7. 重ね合わされた残基の数

1つのタンパク質ペアを例に、このデータ構造に慣れ親しんでみます。

print("Structures in the data set:")
cols = align_df_proteins.columns
rows = align_df_proteins.index
print(cols)
print(rows)

example_col = cols[0]
example_row = rows[2]
print("\nExample align return values for the {}-{} pair: ".format(example_col, example_row))
print(align_df_proteins.loc[example_col, example_row])

print("\nRMSD value (in Angström) after refinement for the {}-{} pair: ".format(example_col, example_row))
print(align_df_proteins.loc[example_col, example_row][0])
"""
    Structures in the data set:
    Index(['1XBB', '2PL0', '4R7I', '3HEC', '4CSV', '1T46', '3FW1'], dtype='object')
    Index(['1XBB', '2PL0', '4R7I', '3HEC', '4CSV', '1T46', '3FW1'], dtype='object')
    
    Example align return values for the 1XBB-4R7I pair: 
    (1.0288203954696655, 1167, 5, 6.308946132659912, 1539, 404.5, 245)
    
    RMSD value (in Angström) after refinement for the 1XBB-4R7I pair: 
    1.0288203954696655
"""

さて、全てのペアについてリファイン後のRMSD(あるいはその他のalignの戻り値のどれでも)を取り出し、さらなる解析のためPandsデータフレームとします。

従って、pymol_align_resultsの全てのペアに関するalign関数の戻り値、および、知りたいpositionの戻り値の位置を入力とする関数を定義します。

def extract_align_info(align_df, position):
    if position in [0, 3, 5]:
        return align_df.applymap(lambda x: np.float(x[position]))
    if position in [1, 2, 4, 6]:
        return align_df.applymap(lambda x: np.int(x[position]))
    else:
        print("Position not available.")

例えば、タンパク質の各ペアについて、重ね合わされた残基の数をチェックすることができます(alignの戻り値の最後の位置)。

extract_align_info(align_df_proteins, 6)
1XBB 2PL0 4R7I 3HEC 4CSV 1T46 3FW1
1XBB 292 272 245 159 249 258 81
2PL0 272 290 267 207 269 261 146
4R7I 245 267 344 186 263 308 26
3HEC 159 207 186 350 243 196 90
4CSV 249 269 263 243 276 269 47
1T46 258 261 308 196 269 316 101
3FW1 81 146 26 90 47 101 238

ここでは、ほとんどのタンパク質のペアについて、残基の大部分が重ね合わせできていることがわかります(ここではEGFR配列の長さは270から350の範囲をとっています)。しかしながら、3FW1の配列のアラインメントは低くなっています。

参考までに:全てのペアの配列アラインメントは../data/T10/alignment/でチェックできます。例を下に示します。

example_alignment_path = glob.glob("../data/T10/alignment/*")[3]
f = open(example_alignment_path)
file_content = f.readlines()
for i in file_content:
    print(i, end="")

"""
  CLUSTAL
    
    i            MALEEIRPKEVYLDRKLLTLED---------------------------------KELGSGNF
    m            MKKGHHHHHHGQKPKYQVRWKIIESYEGNSYTFIDPTQLPYNEKWEFPRNNLQFGKTLGAGAF
                                                                        *.**.*.*
    
    i            GTVKKGYYQM----KKVVKTVAVKILKNEANDPALKDELLAEANVMQQL-DNPYIVRMIGICE
    m            GKVVEATAFGLGKEDAVLK-VAVKMLKSTAHADE-KEALMSELKIMSHLGQHENIVNLLGACT
                 *.*..         .     ****.**        *..*..*...*..* ....**...*.* 
    
    i            AESWMLVM-EMAELGPLNKYLQQNR-----------------------HVKDKNIIELVHQVS
    m            HGGPVLVITEYCTYGDLLNFLRRKAEAMLGPSLSPGQDPEGLDKEDGRPLELRDLLHFSSQVA
                          *....*.*...*...                        ............**.
    
    i            MGMKYLEESNFVHRDLAARNVLLVTQHYAKISDFGLSKALRADENYYKAQTHGKWPVKWYAPE
    m            QGMAFLASKNCIHRDVAARNVLLTNGHVAKIGDFGLARDIMNDSNYIVKGNA-RLPVKWMAPE
                 .**..*...*..***.*******.   .***.*                    ..****.***
    
    i            CINYYKFSSKSDVWSFGVLMWEAFSYGQKPYRGM-KGSEVTAMLEKGERMGCPAGCPREMYDL
    m            SIFDSVYTVQSDVWSYGILLWEIFSLGLNPYPGILVNSKFYKLVKDGYQMAQPAFAPKNIYSI
                 .*........*****.*.*.**.**.*..**.*  .   .......*..*..**..*...*..
    
    i            MNLCWTYDVENRPGFAAVELRLRNYYYDVVNEGHHHHHH                        
    m            MQACWALEPTHRPTFQQITSFLQEQAQEDRR--------                        
                 *..**......**.*...                                             
"""

次のステップでは、リファイン後のRMSD値のデータフレームを作成します(alignの戻り値の最初の位置)。

rmsd = extract_align_info(align_df_proteins, 0)
rmsd
1XBB 2PL0 4R7I 3HEC 4CSV 1T46 3FW1
1XBB 0.000000 1.549436 1.028820 2.560788 0.897684 1.132188 15.648476
2PL0 1.549436 0.000000 1.161164 1.894447 0.871177 1.021322 19.819407
4R7I 1.028820 1.161164 0.000000 2.022787 0.736076 0.687878 7.503986
3HEC 2.560788 1.894447 2.022787 0.000000 1.683216 1.912085 14.695195
4CSV 0.897684 0.871177 0.736076 1.683216 0.000000 0.671700 16.296894
1T46 1.132188 1.021322 0.687878 1.912085 0.671700 0.000000 15.937148
3FW1 15.648476 19.819407 7.503986 14.695195 16.296894 15.937148 0.000000

この改善されたRMSDの結果をヒートマップとして可視化します。

sns.heatmap(rmsd, linewidths=1, annot=True, cbar_kws={"label": "RMSD ($\AA$)"}, cmap="Blues")
plt.show()

f:id:magattaca:20200506010049p:plain

リファイン後のRMSDに基づいてタンパク質の類似性を見るためにヒートマップをクラスタリングします。

def plot_clustermap(rmsd, title):
    g = sns.clustermap(rmsd, linewidths=1, annot=True, cbar_kws={"label": "RMSD ($\AA$)"}, cmap="Blues")
    plt.setp(g.ax_heatmap.get_xticklabels(), rotation=0)
    plt.setp(g.ax_heatmap.get_yticklabels(), rotation=0)
    sns.set(font_scale=1.5)
    
    # プロットの保存ーテキストボックスを含めるためにbbox_inchesを使います:
    # https://stackoverflow.com/questions/44642082/text-or-legend-cut-from-matplotlib-figure-on-savefig?rq=1
    plt.savefig("../data/T10/bsc_{}.png".format(title), dpi=300, bbox_inches="tight", transparent=True)

    plt.show()
plot_clustermap(rmsd, "protein")

f:id:magattaca:20200506010140p:plain

RMSDの比較により、一つのタンパク質、つまり3FW1が他の全てのタンパク質とは異なっていることが示されました(既にアラインメントを3次元で目視で検証した結果と、重ね合わされた残基の数に基づいておこなった議論のとおりです)。

タンパク質は触媒する化学反応によって、いわゆるEC(Enzyme Commission)番号で分類されます。ここでは、このEC番号をタンパク質が属する酵素グループを確認するために使います。

# PDBからPDB IDのEC番号を取得
pdb_all_info = [get_all_info(pdb_id) for pdb_id in pdb_ids]
ec_numbers = [i["polymer"]["enzClass"]["@ec"] for i in pdb_all_info]
target_set = {"pdb_id": pdb_ids,
              "ec_number": ec_numbers}
target_set = pd.DataFrame(target_set)
target_set
pdb_id ec_number
0 2PL0 2.7.10.2
1 3FW1 1.10.5.1
2 3HEC 2.7.11.24
3 4R7I 2.7.10.1
4 1XBB 2.7.10.2
5 4CSV 2.7.10.2
6 1T46 2.7.10.1

3FW1、つまりヒトキノン還元酵素2(human quinone reductase 2、NQO2)はEC class 1、つまり酸化還元酵素(oxidoreductases)に属している一方で、他のタンパク質は全てEC class 2.7、つまりリン酸転移酵素phosphorus transferases)に属しています。EC class 2.7はチロシンキナーゼ(EC 2.7.10.2)を含んでおり、これはイマチニブの目標とするターゲットです。3FW1は「医薬品デザインと慢性骨髄性白血病(chronic myelogenous leukemia)患者の治療に潜在的に関与しうる」オフターゲットとして報告がなされています(BMC Struct. Biol. (2009), 9, 10.1186/1472-6807-9-7)。

PDBの構造をアラインメント(結合サイト)

ここまでアラインメントとRMSDのリファインメントにタンパク質全体の構造を使ってきました。ですが、リガンドはタンパク質の結合サイトだけに結合し、従ってタンパク質全体よりもむしろ結合サイトの類似性がオフターゲットを予測するうえでより確からしい基準とされています。

リガンド原子のいずれかの10 Å 以内となる残基全てを選択することで、タンパク質の結合サイトを定義します。これらの結合サイトの残基をアラインメントに用い、 Cɑ 原子(タンパク質主鎖)だけをRMSDのリファインメントに用います。ここでは、pdb_idsリストの全ての構造を一番最初の構造に対してアラインメントをとります。

# PyMolの再初期化
objPMV.server.do("reinitialize")

# タンパク質ファイルの読み込み
for pdb_id in pdb_ids:
    cmd = "fetch " + pdb_id
    objPMV.server.do(cmd)

# cartoonとしてタンパク質を表示(使用しているpymolのversionに応じて必要な可能性があります)
#objPMV.server.do('cmd.show("cartoon","all")')

# オブジェクトを隠す
objPMV.server.do("hide polar_contacts")

# 背景を白色に設定
objPMV.server.do("bg_color white")
# 水分子とイオンを除去
objPMV.server.do("remove solvent")

# 結合サイトのアラインメント
immobile_pdb_id = pdb_ids[0]
for mobile_pdb_id in pdb_ids[1:]:
    # STIのいずれかの原子と一定の半径以内にある原子を選択し、選択を残基全体に拡張
    objPMV.server.do("select mobile_bs, byres " + mobile_pdb_id + " within 10 of (" + mobile_pdb_id + " and resn STI)")
    objPMV.server.do("select immobile_bs, byres " + immobile_pdb_id + " within 10 of (" + immobile_pdb_id + " and resn STI)")
    # アラインメントを実行
    objPMV.server.do("align mobile_bs, immobile_bs")

# 中心に置き、ズーム
objPMV.server.do("center all")
objPMV.server.do("zoom all")
objPMV.server.do("ray 400,400")

# Jupyter notebookにPyMolフレームを表示
objPMV.GetPNG(h=500)

f:id:magattaca:20200506010215p:plain

ペアワイズRMSDの取得(結合サイト)

タンパク質構造全体のアラインメントとRMSDのリファインメントで先に示したのと同様に、ここでは全てのタンパク質の結合サイトのペアについてアラインメントとRMSDを計算するPyMol Pythonスクリプトを実行します。

使用するPyMolコマンドはPyMol Pythonスクリプトの中で説明されています。

f = open("./pymol_scripts/pymol_align_bindingsites.py")
file_content = f.readlines()
for i in file_content:
    print(i, end="")

"""
    '''
    This script aligns all pairs of input proteins (binding sites).
    '''
    # pymolのpathを追加するために以下2行追記
    import sys
    sys.path.append('/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages')
    # 追記終わり
    import pymol
    import pickle
    import pandas as pd
    import sys
    import numpy as np
    
    # Launch PyMol
    pymol.finish_launching()
    
    # Get the parameter (radius from ligand)
    radius = sys.argv[1]
    
    # Load PDB IDs
    pdb_ids = pickle.load(open("../data/T10/pdb_ids.p", "rb"))
    
    # Create DataFrame to store align results
    align_df = pd.DataFrame(index=pdb_ids, columns=pdb_ids)
    
    def align_structures(m1, m2, align_df):
      # Fetch PDB structures
      pymol.cmd.fetch(m1, "i")
      pymol.cmd.fetch(m2, "m")
    
      # Select binding site
      pymol.cmd.select("i_bs", "byres i within " + radius + " of (i and resn STI)")
      pymol.cmd.select("m_bs", "byres m within " + radius + " of (m and resn STI)")
    
      # Align binding sites
      aln = pymol.cmd.align(mobile="m_bs", target="i_bs", quiet=1)
    
      # Save align results to DataFrame
      align_df.loc[m1, m2] = aln
    
      # Reinitialize PyMol for next structure pair to be loaded
      pymol.cmd.reinitialize()
    
      return align_df
    
    # Iterate over all structure pairs
    for n, immobile in enumerate(pdb_ids):
      for mobile in pdb_ids[n:]:
          align_df = align_structures(mobile, immobile, align_df)
          align_df = align_structures(immobile, mobile, align_df)
    
    # Quit PyMol
    pymol.cmd.quit()
    
    # Save results to file
    pickle.dump(align_df, open("../data/T10/align_df_bindingsites.p", "wb"))
"""

PyMol Pythonスクリプトをターミナル経由で実行します(結合サイトの原子を定義するためリガンドータンパク質半径を入力とます)。

# PyMolのalign関数を追加って結合サイトの比較を実施
os.popen("python ./pymol_scripts/pymol_align_bindingsites.py 10")

結合サイト比較のためにalignデータフレームを読み込みます。

align_df_bindingsites = pickle.load(open("../data/T10/align_df_bindingsites.p", "rb"))

データフレームからRMSD値を抜き出します。

rmsd_bs = extract_align_info(align_df_bindingsites, 0)
extract_align_info(align_df_bindingsites, 6)
2PL0 3FW1 3HEC 4R7I 1XBB 4CSV 1T46
2PL0 87 6 67 81 51 79 85
3FW1 6 38 4 6 7 8 6
3HEC 67 4 83 68 48 67 69
4R7I 81 6 68 90 59 83 86
1XBB 51 7 48 59 65 57 55
4CSV 79 8 67 83 57 83 79
1T46 85 6 69 86 55 79 98

RMSDの結果のクラスタリングしたヒートマップを表示します。

# ペアワイズRMSD値をクラスタリングしたヒートマップとして表示
plot_clustermap(rmsd_bs, "bs")

f:id:magattaca:20200506010314p:plain

アラインメントをとった結合サイトのRMSD値はデータセット(EC number 2.7)における3FW1(EC number 1.10.5.1)の非類似性を示していて、PyMol上の目視による検証でもSTIがタンパク質表面上に結合していることがわかります。1XBB-4CSVと1XBB-3HECのペアも類似性が低いことを示している一方で、データセットの残りはRMSD値が低くなっています。

ここで計算したRMSD値は残基の選択(結合サイトの定義)と、あらかじめ与えられた配列のアラインメントの質に依存します。

# ディレクトリを綺麗にする(ダウンロードしたPDBとPyMolアラインメントファイルを除去)
os.popen("rm ./*.cif")
os.popen("rm ../data/T10/alignment/*.aln")

ディスカッション

このトークトリアルでは、イマチニブに結合するタンパク質のセットにおける類似性と非類似性を評価するために、タンパク質全体と結合サイトについて、配列のアラインメントとそれに続くRMSDのリファインメントを行いました。しかしながら、イマチニブのオフターゲットの予測のためには、解かれた構造の巨大なデータベース(PDB)と、イマチニブの目的のターゲット(チロシンキナーゼ)の結合サイトとを比較する必要があります。この場合、類似性の低い配列の比較もおこなうことになるので、より洗練された検索を可能とするために、配列に依存しないアラインメントアルゴリズムを使う手法や、結合サイトの物理化学的特性を含む手法といった、より洗練された手法が求められます。

クイズ

  1. 医薬品のオンターゲット(on-target)とオフターゲット(off-target)という用語を説明してください。
  2. なぜ結合サイト類似性が、クエリタンパク質のオフターゲットを見つけるのに使えるか説明してください?
  3. オフターゲットを予測する上で(i)タンパク質全体と(ii)タンパク質結合サイトのRMSD値がどう役にたつか議論してください。
  4. 結合サイトの情報にについて代わりとなるアプローチについて考えてください(結合サイトの比較のためにどうやって結合サイトをエンコードするか?)。

訳以上

追記

今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。

今回はオフターゲットの予測を目的としたタンパク質の構造比較でした。内容としては計算リソースの課題もあり、同一リガンドに結合することが確からしい(共結晶構造が取得されるレベルの情報がある)タンパク質の比較(RMSD)でした。検証したエントリーの中にEC番号がキナーゼですらないタンパク質も含まれており、こんなタンパク質も共結晶が取られているんだ!ということにびっくりしました。これを実際にPDBの構造全体に行おうとしたらとんでもなく大変なのではないでしょうか??

TeachOpenCADDを開発されたVolkamer研究室は研究テーマとして「オフターゲット予測」をあげていらっしゃいますので、おそらく効率的な検索手法の開発をされているのだと思います。不勉強で研究内容の中身まで踏み込めておりません。すみません。

昨今では配列からのタンパク質構造の予測についてDeep Learningを用いた手法(AlphaFold )の精度向上も話題になっています。配列の1Dアラインメントのみから2D、3Dの構造予測を踏まえてオフターゲットを予測する、なんてこともできるようになるのかもしれませんね(?)。