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の取得(結合サイト)
レファレンス
結合サイト比較:
- 結合サイト比較についてのレビュー: (Curr. Comput. Aided Drug Des. (2008), 4, 209-20) および (J. Med. Chem. (2016), 9, 4121-51)
- PyMolのコマンド
align
についてのドキュメンテーション (PyMolWiki:align
) - Wikipedia記事:平均二乗偏差(RMSD)(Wikipedia: RMSD)および構造の重ね合わせ(superposition)(Wikipedia: structural superposition)
- 構造の重ね合わせ (Book chapter: Algorithms, Applications, and Challenges of Protein Structure Alignment in Advances in Protein Chemistry and Structural Biology (2014), 94, 121-75)
イマチニブ(Imatinib):
- イマチニブについてのレビュー (Nat. Rev. Clin. Oncol. (2016), 13, 431-46)
- イマチニブの選択性の低さ(Promiscuity of imatinib) (J. Biol. (2009), 8, 10.1186/jbiol134)
- Imatinibに関するChEMBLの情報 (ChEMBL: Imatinib)
- ImatinibのPDB上の情報 (PDB: STI)
- イマチニブの副作用 (Drugs.com: Imatinib)
- イマチニブの副作用 (BMC Struct. Biol. (2009), 9, 10.1186/1472-6807-9-7)
理論
オフターゲットタンパク質
オフターゲットは、医薬品あるいは医薬品の代謝物(の一つ)と相互作用するタンパク質で、ターゲットタンパク質として意図されていないものです。オフターゲットによって起こされる分子の応答は、望まぬ副作用につながる可能性があり、あまり害の無いものから非常に有害な副作用まであります。オフターゲットが生じる主な理由は、本来のターゲット(on-target)とオフターゲット(off-target)が互いの結合サイトに類似の構造的モチーフを共に有しており、それゆえ類似のリガンドに結合することができる、ということです。
コンピューターによるオフターゲット予測:結合サイト比較
コンピューター支援による潜在的なオフターゲットの予測の目的は、潜在的な危険性をもつ医療用物質を開発するリスクを最小にすることです。結合サイトの類似性評価にはいくつかのアルゴリズム的アプローチがありますが、常に次の3つの主要なステップで構成されています。
- 結合サイトの符号化(エンコーディング):結合サイトを様々な記述子テクニックでエンコードし、ターゲットデータベースに格納する
- 結合サイトの比較:様々な類似性指標を使って、クエリの結合サイトをターゲットデータベースと比較する
- ターゲットランキング:適切なスコア化に基づきターゲットを順位づけする
様々な類似性指標と既存のツールについてのより詳細な情報は、結合サイト比較に関する次の素晴らしいレビュー2本を参照して下さい(Curr. Comput. Aided Drug Des. (2008), 4, 209-20 そしてJ. Med. Chem. (2016), 9, 4121-51)。
類似性の簡単な指標としてのペアワイズ平均二乗偏差
類似性のスコア化のための簡単かつ直接的な手法は平均二乗偏差(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)
STIの3次元構造を調べるために、オープンソースツールのPyMolを使います(トークトリアルT8のイントロダクションを参照してください)。PyMolでSTIを眺める前に、3次元座標を計算する必要があります。
まず、化合物に水素原子を追加しますが、水素原子はSMILES形式でいつも明示的に記述されているとは限りません。次に、距離幾何学(ディスタンスジオメトリー法、Distance Geometry)を使って化合物の初期座標を取得し、UFF(Universal Force Field)力場を使って化合物構造の最適化を行います。
sti_mol = Chem.AddHs(sti) sti_mol
AllChem.EmbedMolecule(sti_mol) AllChem.UFFOptimizeMolecule(sti_mol) sti_mol
これで、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)
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データセットをフィルタリング
以下のクライテリアに従ってデータセットをフィルタリングするために、pypdb
のget_entity_info
関数を使ってPDB構造のメタ情報を取得します。
- 実験手法でフィルタリング (
xray
) - 分解能でフィルタリング (3 Å以下)
- (簡単にするため)1本鎖のPDB構造のみを残す(single chain)
- Imatinibに結合した構造のみを残す(例、Imatinibに関連づけられているが、結合はしていないPDB構造が除かれる)
- 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.df
DataFrameの辞書となっており、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)
このイメージは美しくカラフルでくるくると巻いていますが、まだ有益な情報を与えてくれません。次のステップではこの構造を互いに重ね合わせます(アラインメント)。
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)
タンパク質のうちの一つ、3FW1は他のタンパク質と比較して重なりが悪くなっています。重なりの良いタンパク質を表示させるために、このタンパク質を隠します。
objPMV.HideObject("3FW1") objPMV.server.do("center all") objPMV.server.do("zoom all") objPMV.server.do("ray 400,400") objPMV.GetPNG(h=500)
多くの部分(例、ヘリックス)では構造の重なりの程度は高いですが、他の部分では重なりは低いか悪くなっています。
ペアワイズ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項目のタプル、を含むデータフレームです:
- リファイン後のRMSD
- リファイン後の重ね合わされた原子の数
- リファインメントサイクルの試行数
- リファイン前のRMSD
- リファイン前の重ね合わされた原子の数
- アラインメントスコアの生の値
- 重ね合わされた残基の数
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()
リファイン後の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")
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)
ペアワイズ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")
アラインメントをとった結合サイトの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)と、イマチニブの目的のターゲット(チロシンキナーゼ)の結合サイトとを比較する必要があります。この場合、類似性の低い配列の比較もおこなうことになるので、より洗練された検索を可能とするために、配列に依存しないアラインメントアルゴリズムを使う手法や、結合サイトの物理化学的特性を含む手法といった、より洗練された手法が求められます。
クイズ
- 医薬品のオンターゲット(on-target)とオフターゲット(off-target)という用語を説明してください。
- なぜ結合サイト類似性が、クエリタンパク質のオフターゲットを見つけるのに使えるか説明してください?
- オフターゲットを予測する上で(i)タンパク質全体と(ii)タンパク質結合サイトのRMSD値がどう役にたつか議論してください。
- 結合サイトの情報にについて代わりとなるアプローチについて考えてください(結合サイトの比較のためにどうやって結合サイトをエンコードするか?)。
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回はオフターゲットの予測を目的としたタンパク質の構造比較でした。内容としては計算リソースの課題もあり、同一リガンドに結合することが確からしい(共結晶構造が取得されるレベルの情報がある)タンパク質の比較(RMSD)でした。検証したエントリーの中にEC番号がキナーゼですらないタンパク質も含まれており、こんなタンパク質も共結晶が取られているんだ!ということにびっくりしました。これを実際にPDBの構造全体に行おうとしたらとんでもなく大変なのではないでしょうか??
TeachOpenCADDを開発されたVolkamer研究室は研究テーマとして「オフターゲット予測」をあげていらっしゃいますので、おそらく効率的な検索手法の開発をされているのだと思います。不勉強で研究内容の中身まで踏み込めておりません。すみません。
昨今では配列からのタンパク質構造の予測についてDeep Learningを用いた手法(AlphaFold )の精度向上も話題になっています。配列の1Dアラインメントのみから2D、3Dの構造予測を踏まえてオフターゲットを予測する、なんてこともできるようになるのかもしれませんね(?)。
TeachOpenCADD トピック9 〜リガンドベースファーマコフォア〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのREADMEをもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!
以下、訳。
- トークトリアル 9
- リガンドベースファーマコフォア(Ligand-based pharmacophores)
- このトークトリアルの目的
- 学習の目標
- レファレンス
- 理論
- 実践
- nglviewで可視化
- ディスカッション
- クイズ
- 追記
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
トークトリアル 9
リガンドベースファーマコフォア(Ligand-based pharmacophores)
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Pratik Dhakal, Florian Gusewski, Jaime Rodríguez-Guerra and Dominique Sydow
注: このノートブックのセルは一つずつ順番に実行してください。全てのセルを一度に実行することも可能ですが、NGLviewのいくつかは意図した通りに動かないかもしれません。
このトークトリアルの目的
このトークトリアルでは、以前のトークトリアルで選択し、重ね合わせを行った既知のEGFRリガンドを用いて、各リガンドのドナー、アクセプター、そして疎水性ファーマコフォアフィーチャーを見つけます。次に、これらのフィーチャーをクラスター化し、組み合わせファーマコフォア(ensemble pharmacophore)を定義します。これは既知のEGFRリガンドのセットの特性を表し、バーチャルクリーニングによる新奇なEGFRリガンドの探索に用いることができます。
学習の目標
理論
- ファーマコフォアモデリング
- 構造ベース(structure-based)とリガンドベース(ligand-based)ファーマコフォアモデリング
- ファーマコフォアを用いたバーチャルスクリーニング
- クラスタリング:k平均法(k-means)
実践
- 先のトークトリアルであらかじめ重ね合わせたリガンドを取得
- NGLViewでリガンドを表示
- ファーマコフォアフィーチャー(pharmacophore feature)を抽出
- 全リガンドのファーマコフォアフィーチャーを表示
- 水素結合ドナー
- 水素結合アクセプター
- 疎水性コンタクト
- フィーチャータイプ毎にフィーチャーの座標を収集
- 組み合わせファーマコフォアの作成
- クラスターの表示
- 水素結合ドナー
- 水素結合アクセプター
- 疎水性コンタクト
- 組み合わせファーマコフォアの表示
レファレンス
- IUPACにおけるファーマコフォアの定義 (Pure & Appl. Chem (1998), 70, 1129-43)
- LigandScoutにおける3Dファーマコフォア (J. Chem. Inf. Model. (2005), 45, 160-9)
- 書籍の章: Pharmacophore Perception and Applications (Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 259-82)
- 書籍の章: Structure-Based Virtual Screening (Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 313-31).
- Monty Kierとファーマコフォアの概念の起源 (Internet Electron. J. Mol. Des. (2007), 6, 271-9)
- Nik StieflによるRDKitを使ったファーマコフォアモデリングのデモンストレーション (RDKit UGM 2016 on GitHub)
理論
ファーマコフォア
コンピューター支援医薬品デザイン(computer-aided drug design)において、ファーマコフォアを用いた医薬品とターゲット分子との相互作用の表現はよく確立された手法です。ファーマコフォアという用語は1998年IUPACの作業部会によって定義されました。即ち、「ファーマコフォアは、特定の生物学的標的の構造との最適な分子間相互作用を確実なものとし、その生物学的応答を引き起こす(あるいは阻害する)ために必要な、立体的、電子的特徴の組み合わせ」です。 (Pure & Appl. Chem. (1998), 70, 1129-43)
言い換えれば、ファーマコフォアはいくつかのファーマコフォアフィーチャーで構成されていて、ファーマコフォアフィーチャーは研究対象のターゲット分子に結合することが観察されているリガンドの重要な立体的、物理化学的特性を表現します。そのような物理化学的特性(フィーチャータイプとも呼ばれます)としては水素結合ドナー/アクセプター、疎水性/芳香族性相互作用、あるいは、正/負電荷を帯びた官能基があり、立体的特性はこれらのフィーチャーの3次元配置によって定義されます。
構造ベースファーマコフォアモデリングとリガンドベースファーマコフォアモデリング
ファーマコフォアモデリングには、生物学的な問いと、手に入るデータソースに応じて、2つの主要なアプローチが使われます。即ち、構造ベースファーマコフォアモデリングとリガンドベースファーマコフォアモデリングです。
構造ベースファーマコフォアモデル(Structure-based pharmacophore models) はタンパク質-リガンド複合体から作成されます。フィーチャーはタンパク質とリガンドの間で観察される相互作用によって定義され、リガンドの結合に関与することが示されているこれらのリガンド部位だけがバーチャルスクリーニングに使われることを保証します。 しかしながら、タンパク質-リガンド複合体の構造は全ての標的タンパク質で手に入る訳ではありません。この場合、例えばドッキングによってリガンドを標的の結合サイトにモデリングすることで複合体の構造を生成するか、あるいは、タンパク質とリガンドが相互作用する可能性がある場所を見つけるために、標的タンパク質の結合サイトだけを使ってファーコフォアモデリングを行います。
リガンドベースファーマコフォアモデル(Ligand-based pharmacophore models) は研究対象の標的分子に結合することが知られている一連のリガンドに基づきます。これらのリガンドに共通してみられる化学的特徴でファーマコフォアモデルを作成します。この手法は複数の既知リガンドがある標的タンパク質で、タンパク質ーリガンド複合体構造が無い場合に使われます。このトークトリアルでは、既知のEGFRリガンドセットを使ってリガンドベースのファーマコフォアモデリングを行います。
ファーマコフォアモデリングに関してさらに知りたい場合は、 (Pharmacophore Perception and Applications: Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 259-82) and (J. Chem. Inf. Model. (2005), 45, 160-9)をお勧めします。
ファーマコフォアを用いたバーチャルスクリーニング
先にトークトリアル 4で説明したように、バーチャルクリーニング(Virtual Screening, VS)は、(クエリによって代表される)研究対象のターゲット分子に結合する可能性が最も高い(ライブラリ中の)低分子を見つけるために、巨大な化合物ライブラリに対してクエリ(例えばこのトークトリアル 9ではファーマコフォアモデル、トークトリアル 4ではクエリ化合物)のスクリーニングを実施することです。ファーマコフォアベースのバーチャルスクリーニングでは、化合物ライブラリがファーマコフォアモデルに対して化合物毎に突き合わせられ、最も整合性の良かった結果にもとづきランク付けされます(Structure-Based Virtual Screening: Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 313-31)。
クラスタリング:k平均法
このトークトリアルでは、いくつかのリガンドベースファーマコフォアのフィーチャーポイントをクラスタリングすることで組み合わせファーマコフォアを作成します。クラスタリングのアルゴリズムとして、データセットをk個のクラスターに分類するために使われる、k平均法クラスタリングを使います:
- k個の異なる重心(centroid)を選択し、データセットの各ポイントを最も近い重心に割り当てる
- 現在のクラスターに基づいて新しい重心を計算し、データセットの各ポイントが新しい重心のうち最も近いものに割り当てられる
- 重心が安定化するまでこの過程を繰り返す
実践
import os, glob # RDKit from rdkit import RDConfig, Chem, Geometry, DistanceGeometry from rdkit.Chem import ChemicalFeatures, rdDistGeom, Draw, rdMolTransforms, AllChem from rdkit.Chem.Draw import IPythonConsole, DrawingOptions from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib from rdkit.Numerics import rdAlignment IPythonConsole.ipython_useSVG=True import collections import pandas as pd import math from sklearn import datasets, cluster import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from collections import Counter # For handling the labels import operator
先のトークトリアルのあらかじめ重ね合わせたリガンドを取得
先のトークトリアルで重ね合わせたリガンドを取得します。
まず、全リガンドPDBファイルへのファイルパスを取得します。
mol_files = list(glob.glob("../data/T8/*_lig.pdb")) mol_files # ['../data/T8/5UG8_lig.pdb', # '../data/T8/3POZ_lig.pdb', # '../data/T8/5UG9_lig.pdb', # '../data/T8/5HG8_lig.pdb']
pdb_ids = [os.path.basename(i).split("_")[0] for i in mol_files] pdb_ids # ['5UG8', '3POZ', '5UG9', '5HG8']
訳注(2020/05)
まず、Python標準ライブラリのglob
モジュールを使うことでワイルドカード*
を使って一連のファイルのパスを取得しています。
次にos.path
モジュールのbasename()
で、パス文字列からファイル名を取得し、split
で分割することでPDB ID
だけにしています。
訳注ここまで
次に、RDKitを使ってこれらのPDBファイルから全てのリガンドを読み込みます。
mols = [] for mol_file in mol_files: mol = Chem.MolFromPDBFile(mol_file, removeHs=False) if mol is None: print(mol_file, 'could not be read') else: Chem.SanitizeMol(mol) print(Chem.MolToSmiles(mol)) mols.append(mol) rangeMols = range(1, len(mols)+1) print('Number of molecules: ', len(mols)) # CCC(O)N[C@@H]1CN(C2NC(NC3CNN(C)C3)C3NCN(C(C)C)C3N2)C[C@H]1F # CC(C)(O)CC(O)NCCN1CCC2NCNC(NC3CCC(OC4CCCC(C(F)(F)F)C4)C(Cl)C3)[C@H]21 # CCC(O)N[C@@H]1CN(C2NC(NC3CN(C)NC3OC)C3NCN(C(C)C)C3N2)C[C@H]1F # CCC(O)NC1CCCC(OC2NC(NC3CNN(C)C3)NC3NCCC32)C1 # Number of molecules: 4
Draw.MolsToGridImage(mols, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])
ここで問題に行き当たりました:PDBファイルからリガンドを読み込むと、RDKitはリガンドに対して例えば芳香環の情報を割り当てません。RDKitの関数AssignBondOrdersFromTemplate
を使うと、リファレンス分子に基づいて分子の結合に情報を割り当てます。リファレンス分子として例えば、我々の場合化合物のSMILESパターンを使います。
もっと情報が知りたい場合は(RDKitディスカッションの「PDBの非タンパク質分子の芳香族が認識されない」)と (RDKitドキュメンデーションのAssignBondOrdersFromTemplate
)をチェックしてください。
# PDBリガンド構造のSMILESを読み込む ligs = pd.read_csv("../data/T8/PDB_top_ligands.csv", sep="\t") # pdb_idsと同じ順番でSMILESを取得 ligs_smiles = [ligs[ligs["@structureId"]==pdb_id]["smiles"].values[0] for pdb_id in pdb_ids] # SMILESからRDKit Molオブジェクトを生成 refmols = [Chem.MolFromSmiles(smiles) for smiles in ligs_smiles] # SMILESパターン(refmols)に基づき化合物(mols)に結合次数を割り当てる mols = [AllChem.AssignBondOrdersFromTemplate(refmol, mol) for refmol, mol in zip(refmols, mols)]
Draw.MolsToGridImage(mols, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])
2次元で化合物を見ることもできます(この例では元々の座標を保持しておくために化合物をコピーしておきます)。
mols_2D = [] for mol in mols: tmp=Chem.Mol(mol) AllChem.Compute2DCoords(tmp) mols_2D.append(tmp) Draw.MolsToGridImage(mols_2D, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])
nglview
で可視化
ImportError
あるいはModuleNotFoundError
の例外がでた場合は、下のセルのタイプをCode
に切り替えて実行してください。
!pip install nglview !nglview install
import nglview as nv
import time def show_ligands(rdkit_mols): v = nv.NGLWidget() for mol in rdkit_mols: c = v.add_component(Chem.MolToPDBBlock(mol), ext="pdb") time.sleep(0.1) c.clear() c.add_ball_and_stick(multipleBond=True) return v
v = show_ligands(mols) v
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
ファーマコフォアフィーチャーの抽出
上で述べたように、このトークトリアルの目的はリガンドセットからリガンドベースのファーマコフォアの組み合わせを作成することです。まず最初に、リガンド毎にファーマコフォアフィーチャーを取り出す必要があります。そこで、フィーチャーファクトリーを(デフォルトのフィーチャーの定義で)読み込みます。
RDKitドキュメンテーションの「chemical features and pharmacophores」も参照してください。.
ffact = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'))
RDKitに実装されているファーマコフォアフィーチャーを見てみましょう。
list(ffact.GetFeatureDefs().keys()) """ ['Donor.SingleAtomDonor', 'Acceptor.SingleAtomAcceptor', 'NegIonizable.AcidicGroup', 'PosIonizable.BasicGroup', 'PosIonizable.PosN', 'PosIonizable.Imidazole', 'PosIonizable.Guanidine', 'ZnBinder.ZnBinder1', 'ZnBinder.ZnBinder2', 'ZnBinder.ZnBinder3', 'ZnBinder.ZnBinder4', 'ZnBinder.ZnBinder5', 'ZnBinder.ZnBinder6', 'Aromatic.Arom4', 'Aromatic.Arom5', 'Aromatic.Arom6', 'Aromatic.Arom7', 'Aromatic.Arom8', 'Hydrophobe.ThreeWayAttach', 'Hydrophobe.ChainTwoWayAttach', 'LumpedHydrophobe.Nitro2', 'LumpedHydrophobe.RH6_6', 'LumpedHydrophobe.RH5_5', 'LumpedHydrophobe.RH4_4', 'LumpedHydrophobe.RH3_3', 'LumpedHydrophobe.tButyl', 'LumpedHydrophobe.iPropyl'] """
一例として、参考例の化合物の全てのフィーチャーを取得してみます。
m1 = mols[0] feats = ffact.GetFeaturesForMol(m1) print('Number of features found:',len(feats)) # Number of features found: 14
フィーチャーのタイプ(RDKitではfamilyと呼ばれています)はGetFamily()
で取得することができます。
feats[0].GetFamily() # 'Donor'
参考例の化合物のフィーチャーのタイプについてその頻度を取得します。
feats_freq = collections.Counter([x.GetFamily() for x in feats]) feats_freq """ Counter({'Donor': 2, 'Acceptor': 6, 'PosIonizable': 1, 'Aromatic': 3, 'Hydrophobe': 1, 'LumpedHydrophobe': 1}) """
訳注(2020/05)
Python標準ライブラリcollections
のCounter
クラスにより、与えられたリストの要素の出現頻度について辞書型オブジェクトが生成されています。keyが要素、値が出現頻度で、ここではフィーチャータイプと頻度の辞書になっています。
訳注ここまで
上記の関数をリガンドセットの全化合物に適用します。化合物ごとのフィーチャータイプの頻度をデータフレームとして表示します。
# 化合物ごとのフィーチャータイプの頻度を取得 mols_feats_freq = [] for i in mols: feats = [x.GetFamily() for x in ffact.GetFeaturesForMol(i)] feats_freq = collections.Counter(feats) mols_feats_freq.append(feats_freq) # データをデータフレームとして表示 p = pd.DataFrame(mols_feats_freq, index=["m"+str(i) for i in range(1, len(mols)+1)]).fillna(0).astype(int) p.transpose()
m1 | m2 | m3 | m4 | |
---|---|---|---|---|
Donor | 2 | 3 | 2 | 4 |
Acceptor | 6 | 5 | 7 | 5 |
PosIonizable | 1 | 0 | 1 | 0 |
Aromatic | 3 | 4 | 3 | 4 |
Hydrophobe | 1 | 3 | 1 | 2 |
LumpedHydrophobe | 1 | 2 | 1 | 1 |
この先、このトークとリアルでは次のフィーチャータイプのみに焦点をあてます。すなわち水素結合アクセプター(acceptors)、水素結合ドナー(donors)、そして疎水性コンタクト(hydrophobics)です。
フィーチャータイプごと、そして化合物ごとに、RDKitのフィーチャーオブジェクトを取得します。
acceptors = [] donors = [] hydrophobics = [] for i in mols: acceptors.append(ffact.GetFeaturesForMol(i, includeOnly='Acceptor')) donors.append(ffact.GetFeaturesForMol(i, includeOnly='Donor')) hydrophobics.append(ffact.GetFeaturesForMol(i, includeOnly='Hydrophobe')) features = {"donors": donors, "acceptors": acceptors, "hydrophobics": hydrophobics}
全てのリガンドのファーマコフォアフィーチャーを表示
ファーマコフォアフィーチャーのタイプは大抵の場合、定義された色で表示されます。例えば、水素結合ドナー、水素結合アクセプターそして、疎水性コンタクトはそれぞれ緑、赤、黄色でよく表されます。
feature_colors = {"donors": (0,0.9,0), # Green "acceptors": (0.9,0,0), # Red "hydrophobics": (1,0.9,0)} # Yellow
def visualize_features(molecules, feature_type, features, color="yellow", sphere_radius=0.5): print("Number of", feature_type, "in all ligands:", sum([len(i) for i in features])) v = show_ligands(molecules) for i, feature in enumerate(features, 1): for feat in feature: loc = list(feat.GetPos()) label = f"{feature_type}_{i}" v.shape.add_sphere(loc, color, sphere_radius, label) return v
考察中のフィーチャータイプのフィーチャーを可視化するためにこの関数を使います。
水素結合ドナー
feature_type = "donors" v = visualize_features(mols, feature_type, features[feature_type], feature_colors[feature_type]) v # Number of donors in all ligands: 11
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
水素結合アクセプター
feature_type = "acceptors" v = visualize_features(mols, feature_type, features[feature_type], feature_colors[feature_type]) v # Number of acceptors in all ligands: 23
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
疎水性コンタクト
feature_type = "hydrophobics" v = visualize_features(mols, feature_type, features[feature_type], feature_colors[feature_type]) v # Number of hydrophobics in all ligands: 7
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
フィーチャータイプ毎にフィーチャーの座標を集める
(フィーチャータイプ毎に)フィーチャーをクラスタリングしたいので、(フィーチャータイプ毎に)フィーチャーのすべての座標を集めます。
features_coord = {"donors": [list(item.GetPos()) for sublist in features["donors"] for item in sublist], "acceptors": [list(item.GetPos()) for sublist in features["acceptors"] for item in sublist], "hydrophobics": [list(item.GetPos()) for sublist in features["hydrophobics"] for item in sublist]}
これで全てのフィーチャーの座標が手に入りました。例えばacceptorフィーチャーの全ての座標は次のようになります。
features_coord["acceptors"] """ [[-13.788, 14.818, -27.097], [-12.118, 16.261, -28.026], [-11.376, 12.959, -29.238], [-16.74, 10.887, -25.841], [-16.005, 17.279, -22.5], [-15.388, 19.42, -27.168], [-11.052, 12.941, -29.208], [-12.995, 19.311, -27.816], [-10.269, 15.094, -29.72], [-8.077, 20.488, -31.958], [-12.686, 20.459, -24.958], [-13.863, 14.66, -27.134], [-12.175, 16.118, -27.994], [-11.411, 12.856, -29.266], [-17.026, 10.809, -26.396], [-15.601, 9.901, -28.086], [-16.044, 17.165, -22.506], [-15.445, 19.229, -27.215], [-13.162, 14.493, -28.185], [-11.603, 12.928, -29.136], [-12.896, 16.874, -28.352], [-16.001, 17.24, -22.411], [-16.523, 11.38, -26.098]] """
組み合わせファーマコフォアの作成
組み合わせファーマコフォアを作成するために、k平均法クラスタリングをつかってフィーチャータイプ毎にフィーチャーをクラスタリングします。
k平均法クラスタリングの静的パラメーターを設定
kq
: このパラメーターで、フィーチャーポイントの数に応じてフィーチャータイプ毎のクラスター数 k
を決定します。つまり、フィーチャータイプ毎に:
k
= number_of_features / kq
となります。
# k quotient (kq) はk平均法でkを決めるために使われます: k = number of feature points / kq # kq は全てのクラスターについてk(フィーチャークラスター)が少なくとも1クラスターで、4-5クラスターよりも大きくならないように選択するべきです。 kq = 7
クラスター選別のための静的パラメーターを設定
min_cluster_size
: リガンドの集合のなかのほとんどの化合物がもっているフィーチャーを含む可能性があるクラスターだけを残します。したがって、この変数は、我々のリガンドの集合の化合物数の75%に設定します。
top_cluster_number
: このパラメーターで、最も大きいクラスターだけを選択します。
# クラスタリングの閾値:数値=閾値のパーセンテージ min_cluster_size = int(len(mols) * 0.75) # 上位のフィーチャーだけを表示 top_cluster_number = 4
k平均法クラスタリングとクラスター選択の関数を定義
k平均法クラスタリングから導いたクラスターの中心を計算する関数を定義します。
def clustering(feature_coord, kd): ''' この関数は入力のフィーチャーの座標のk平均法クラスタリングを計算します ''' # パラメーターkを、"k quotient"でフィーチャーの数を割った値として定義 k = math.ceil(len(feature_coord) / kq) k = 2 if k == 1 else k # このトークトリアルの例のにあうように疎水性コンタクトのkを調整 print('Clustering: \nVariable k in k-means: %d of %d points\n'%(k, len(feature_coord))) # k平均法を初期化 k_means = cluster.KMeans(n_clusters=k) # k平均法クラスタリングを計算 k_means.fit(feature_coord) # クラスターを返す return k_means
クラスターをサイズで並べ替え、最も大きいクラスターのインデックスのリストを出力する関数を定義します。
def get_clusters(k_means, min_cluster_size, top_cluster_number): ''' この関数は入力されたk平均法クラスタリングの情報を取得します: * 各フィーチャーのクラスターのラベルを取得 * クラスターのサイズを数え、クラスターサイズによってクラスターのインデックスを並べ替える * サイズに基づきクラスターを選択 * 選択したクラスターのインデックスを返す ''' # サイズによってクラスターを並べ替え、最大のもののみを表示 feature_labels = k_means.labels_ print('Cluster labels for all features: \n%s\n'% feature_labels) feature_labels_count = Counter(feature_labels) print('Cluster label counter: \n%s\n'% feature_labels_count) feature_labels_count = sorted(feature_labels_count.items(), key=operator.itemgetter(1), reverse=True) print('Sorted cluster label counters: \n%s\n'% feature_labels_count) # 閾値よりも大きなクラスターで、最大のものの番号を取得(選択されたクラスター) cluster_indices_sel = [] for cluster_index, cluster_size in feature_labels_count: # feature_labels_count = list of (cluster_index, cluster_size) if cluster_size >= min_cluster_size and top_cluster_number > 0: cluster_indices_sel.append(cluster_index) top_cluster_number -= 1 print('Cluster indices of selected clusters: \n%s\n'% cluster_indices_sel) return cluster_indices_sel
訳注(2020/05)
feature_labels_count
は、collections
のCounter
でクラスターのラベルを集計(要素の出現頻度の辞書を作成)した後に並べ替えています。counter
クラスは辞書型で順番が保証されないので、items()
とsorted()
を使うことでソートされたタプル(key, value)
のリストを得ています。辞書のvalueでソートしたいので、sorted()の引数のkey
に、operator.itemgetter()
を使用しています。(itemgetter(0)
ではなく)itemgetter(1)
とすることで辞書のvalueが使われることになります。
訳注ここまで
フィーチャーをクラスタリング
各フィーチャータイプについて、定義したclustering
関数を使ってk平均法クラスタリングを実施します。
k_means = {"donors": clustering(features_coord["donors"], kq), "acceptors": clustering(features_coord["acceptors"], kq), "hydrophobics": clustering(features_coord["hydrophobics"], kq)} """ Clustering: Variable k in k-means: 2 of 11 points Clustering: Variable k in k-means: 4 of 23 points Clustering: Variable k in k-means: 2 of 7 points """
適切なクラスターの選択
各フィーチャータイプについて、定義したget_clusters
関数をつかって適切なクラスターを選択します。
print("Hydrogen bond donors\n") cluster_indices_sel_don = get_clusters(k_means["donors"], min_cluster_size, top_cluster_number) """ Hydrogen bond donors Cluster labels for all features: [1 0 0 0 1 1 0 1 1 0 1] Cluster label counter: Counter({1: 6, 0: 5}) Sorted cluster label counters: [(1, 6), (0, 5)] Cluster indices of selected clusters: [1, 0] """
print("Hydrogen bond acceptors\n") cluster_indices_sel_acc = get_clusters(k_means["acceptors"], min_cluster_size, top_cluster_number) """ Hydrogen bond acceptors Cluster labels for all features: [0 0 0 2 1 1 0 1 0 3 1 0 0 0 2 2 1 1 0 0 0 1 2] Cluster label counter: Counter({0: 11, 1: 7, 2: 4, 3: 1}) Sorted cluster label counters: [(0, 11), (1, 7), (2, 4), (3, 1)] Cluster indices of selected clusters: [0, 1, 2] """
print("Hydrophobic contacts\n") cluster_indices_sel_h = get_clusters(k_means["hydrophobics"], min_cluster_size, top_cluster_number) """ Hydrophobic contacts Cluster labels for all features: [1 0 0 1 1 0 1] Cluster label counter: Counter({1: 4, 0: 3}) Sorted cluster label counters: [(1, 4), (0, 3)] Cluster indices of selected clusters: [1, 0] """
cluster_indices_sel = {"donors": cluster_indices_sel_don, "acceptors": cluster_indices_sel_acc, "hydrophobics": cluster_indices_sel_h}
選択したクラスターの座標を取得
def get_selected_cluster_center_coords(k_means, cluster_indices_sel, feature_type): ''' この関数は選択したクラスターのクラスター中心の座標を(インデックスで)取得します。 ''' # ある特定のフィーチャータイプについてクラスター中心を取得 cluster_centers = k_means[feature_type].cluster_centers_ # (インデックスにより要素を選択するために)リストに型変換したのち、PandasのSeriesに型変換する cluster_centers = pd.Series(cluster_centers.tolist()) # 選択したクラスターのインデックスでクラスター中心を選択 cluster_centers_sel = cluster_centers[cluster_indices_sel[feature_type]] # リストに型変換しリストを返す return list(cluster_centers_sel)
cluster_centers_sel = {"donors": get_selected_cluster_center_coords(k_means, cluster_indices_sel, "donors"), "acceptors": get_selected_cluster_center_coords(k_means, cluster_indices_sel, "acceptors"), "hydrophobics": get_selected_cluster_center_coords(k_means, cluster_indices_sel, "hydrophobics")}
cluster_centers_sel["acceptors"] """ [[-12.155727272727272, 14.545636363636364, -28.48690909090909], [-14.937714285714286, 18.586142857142857, -24.93914285714286], [-16.4725, 10.744250000000001, -26.605249999999998]] """
クラスターを表示
フィーチャータイプ毎に、全ての化合物と全てのフィチャーポイントと一緒にクラスター中心を可視化します。
def visualize_clusters(molecules, feature_type, features, clusters, color="yellow", feature_radius=0.5, cluster_radius=1): v = visualize_features(molecules, feature_type, features, color=color, sphere_radius=feature_radius) for i, center in enumerate(clusters, 1): v.shape.add_sphere(list(center), color, cluster_radius, f"cluster_{feature_type}_{i}") return v
水素結合ドナー
feature_type = "donors" v = visualize_clusters(mols, feature_type, features[feature_type], cluster_centers_sel[feature_type], feature_colors[feature_type]) v # Number of donors in all ligands: 11
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
水素結合アクセプター
feature_type = "acceptors" v = visualize_clusters(mols, feature_type, features[feature_type], cluster_centers_sel[feature_type], feature_colors[feature_type]) v # Number of acceptors in all ligands: 23
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
疎水性コンタクト
feature_type = "hydrophobics" v = visualize_clusters(mols, feature_type, features[feature_type], cluster_centers_sel[feature_type], feature_colors[feature_type]) v # Number of hydrophobics in all ligands: 7
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
組み合わせファーマコフォアの表示
この最後のステップでは、クラスタリングしたファーマコフォアのフィーチャー(即ち、水素結合ドナー、アクセプターと疎水性コンタクト)を組み合わせて、一つの組み合わせファーマコフォアを作成します。この組み合わせファーマコフォアは4つの選択したリガンドのファーマコフォアの特徴を表します。
v = show_ligands(mols) # クラスターの読み込み for feature_type in cluster_indices_sel.keys(): centers = cluster_centers_sel[feature_type] for i, loc in enumerate(centers): sphere_radius = 1 feature_color = feature_colors[feature_type] label = f"{feature_type}_c{i}" v.shape.add_sphere(loc, feature_color, sphere_radius, label) v
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
ディスカッション
このトークトリアルでは、EGFRに結合することが知られているリガンドのセットを、予め重ね合わせたものを用いて、組み合わせファーマコフォアモデルを作成しました。これで、既知のリガンドに観測された立体的、物理化学的特徴を示し、EGFR結合サイトに結合するかもしれないような新奇な低分子を見つけることを目的として、低分子化合物の巨大なライブラリのバーチャルスクリーニングを行うために、この組み合わせファーマコフォアモデルを使うことができます。
スクリーニングの前に、通常、ファーアマコフォアモデルはさらに最適化されます。例えば、スクリーニングに用いるフィーチャーの特徴の数を減らすために、生物学的な知識(ある相互作用は重要である一方、他の相互作用はそうではないという報告されているかもしれません)、あるいは化学的専門知識に基づいて、いくつかのフィーチャーは除外されるかもしれません。
このトークトリアルではバーチャルスクリーニングは扱いませんが、RDKitを使ってファーマコフォアモデリングとバーチャルスクリーニングをデモンストレーションしている、Nik Stieflによる素晴らしいチュートリアルに言及しておきます。(RDKit UGM 2016 on GitHub).
ファーマコフォアフィーチャーのクラスタリングを行うためにk平均法クラスタリングを用いました。このクラスタリング手法は、使用者が予めクラスターの数を定義する必要があるという欠点があります。通常、クラスターの数は、クラスタリングの実施前(あるいはクラスターの精度を高めていく過程)に、ポイントの分布を目視で調査して決めるので、ファーマコフォアを自動的に作成するには妨げとなります。この解決策として、密度に基づくクラスタリング手法(とk平均法クラスタリングの組み合わせも)が挙げられます。
クイズ
- ファーマコフォアフィーチャーとファーマコフォアを説明してください。
- 構造ベースのファーマコフォアモデリングとリガンドベースのファーマコフォアモデリングの違いを説明してください。
- 組み合わせファーマコフォアを導く方法を説明してください。
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回のトピックに関連して日本語で読める記事として株式会社理論創薬研究所 吉森篤史氏による以下の記事をお勧めします。 RDKitを利用したファーマコフォアの定義とその使用方法について解説がされていてとてもわかりやすいです。
- SAR News No.19 (Oct. 2010) 「Chemoinformatics Toolkits を用いた創薬システム開発における ラピッドプロトタイピング」
欲を言えば作ったファーマコフォアモデルによるバーチャルスクリーニングの実施まで行なってほしかったですが、計算リソースとデータセットの準備だけでも大変そうですね。
勝手にいろいろと引用させていただいておりますので問題があれば削除します。
TeachOpenCADD トピック8 〜タンパク質データの取得〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル8をもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下、訳。
トークトリアル 8
タンパク質データの取得: Protein Data Bank (PDB)
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Anja Georgi, Majid Vafadar and Dominique Sydow
このトークトリアルの目的
このトークトリアルでは、次のトークトリアルのための土台を準備します。次のトークトリアルではEGFRに対するリガンドベースの組み合わせファーマコフォア(ligand-based ensemble pharmacophore)を作成します。したがって準備として(i) PDBデータベースからEGFRの全てのPDB IDを取得し、(ii) X線結晶構造解析による構造で、最も質の良い5つのタンパク質ーリガンド複合体構造を取得、そして、(iii) 5つの構造を互いに3Dで重ね合わせ(アラインメント)、次のトークトリアルで使用するためにリガンドを抜き出し保存します。
学習の目標
理論
実践
- 検索するクエリタンパク質の選択
- クエリタンパク質のPDBエントリーに関する統計値を取得
- クエリタンパク質の全PDB IDを取得
- PDBエントリーのメタ情報を取得
- PDBエントリーのメタ情報をフィルタリングし並べ替え
- 上位の構造からリガンドのメタ情報を取得
- 上位リガンド化合物を描画
- タンパク質ーリガンド IDのペアを作成
- PDB構造ファイルの取得
- PDB構造の重ね合わせ
レファレンス
- Protein Data Bank (PDB website)
- PyPDB pythonパッケージ (Bioinformatics (2016), 32, 159-60)
- PyPDB pythonパッケージドキュメンテーション (PyPDB website)
- PyMol selection algebra (PyMolWiki: selection algebra)
理論
Protein Data Bank (PDB)
Protein Data Bank (PDB) は構造生物学に関する情報のデータベースの中で最も包括的なものの一つで、構造ゲノム科学と医薬品デザインといった構造生物学の分野で重要な情報源です。(PDB Webサイト)
構造のデータは、X線結晶構造解析(最も良く使われている手法)、核磁気共鳴(NMR)そしてクライオ電子顕微鏡(cryo-EM)といった構造決定手法によって取得されています。各エントリーについて、データベースに含まれている情報には(i) タンパク質、リガンド、補因子、水分子、そしてイオンといった原子と、原子同士をつなぐ結合の3次元座標、そして(ii) PDB ID、著者、登録日、用いられた構造決定手法や構造の分解能といった構造情報に関するメタ情報があります。
構造の分解能(resolution)は、集められたデータの質の指標で、単位はÅ(オングストローム)です。この値が小さいほど、構造の質はより良いものとなります。
PDB Webサイトはタンパク質の構造(と、手に入る場合はリガンドとの相互作用)の3次元描画と、構造の質の指標(メトリクス)とを提供しています。例として、上皮増殖因子(EGFR)のPDBエントリー PDB ID 3UG5 で見ることができます。
PyPDB
PyPDBはPDBのためのPythonプログラミングインターフェースで、Python3でのみ機能します。 このパッケージを利用することで、バイオインフォマティクスのワークフローにPDB自動検索を組み込み、既存の検索結果に基づいて多数の検索をおこなうプロセスを簡略化することができます。 また、PDBエントリーの情報について高度な検索を行うことも可能になります。PDBは現在RESTful APIを使用しており、標準的なHTMLボキャブラリーを使用して情報を取得することができます。PyPDBはこれらのオブジェクトをXML文字列に変換します。 (Bioinformatics (2016), 32, 159-60)
提供されている機能(関数)のリストはPyPDBドキュメンテーションのWebサイトで見ることができます(PyPDB website)。
実践
# 必要なライブラリのインポート from pypdb import * from pymol import * from rdkit import Chem from rdkit.Chem.Draw import IPythonConsole from rdkit.Chem import Draw from rdkit.Chem import PandasTools IPythonConsole.ipython_useSVG=True import pprint import glob import pandas as pd from array import array import numpy as np import collections import matplotlib.pyplot as plt %matplotlib inline
検索するクエリタンパク質の選択
このトークトリアルではEGFRをクエリタンパク質として使います。EGFRのUniProt IDはP00533
で、以降のPDBデータベース検索でこのIDを使います。
クエリタンパク質のPDBエントリーに関する統計値を取得
まず最初の質問です。「EGFRについて、毎年いくつのエントリーがPDBに登録されていて、全体ではいくつになるのでしょうか?」
PDB Webサイト で検索ターム「P00533
」で検索することができます。2018年10月現在、PDBでの検索結果は179でした。
訳注(2020/05)
2020年5月現在、検索結果は207でした。
訳注ここまで
pypdb
を使うと、PDBデータベースから全てのEGFR構造の登録日(deposition date)を取得することができます。find_dates
関数のパラメーターmax_results
を設定するために、登録されている構造の数が必要となります。
# 注: max_results パラメーターのデフォルトは100でEGFRのエントリーに対して少なすぎます。 # max_results > EGFR構造の最大の数 となるとエラーが出ます。 # 従って実行前に結果がいくつあるかチェックしました。(#179) # 訳注: 2020年5月時点では207で、以下のコードも合わせて変更しています。207ではStopIterationとなったので205としています。 # このデータベースクエリには時間がかかるかもしれません(数分程度) # all_dates = find_dates("P00533", max_results=179) all_dates = find_dates("P00533", max_results=205)
print("Number of EGFR structures found: " + str(len(all_dates))) # Number of EGFR structures found: 205
# 登録日の例 all_dates[:3] # ['2002-03-28', '2002-06-17', '2002-06-17']
登録日から登録年の情報を抽出し、各年毎の登録数のヒストグラムを計算します。
# 登録年の抽出 all_dates = np.asarray(all_dates) all_years = np.asarray([int(depdate[:4]) for depdate in all_dates]) # ヒストグラムの計算 bins = max(all_years)-min(all_years) # ビンの数 = 登録年の範囲 subs_v_time = np.histogram(all_years, bins) # 全エントリー(2018年を除く)のプロット # 訳注: 最新年を除く設定なので用いたデータにより除かれる年は異なります dates, num_entries = subs_v_time[1][:-1], subs_v_time[0] # ヒストグラムの表示 fig = plt.figure() ax = plt.subplot(111) ax.fill_between(dates, 0, num_entries) ax.set_ylabel("New entries per year") ax.set_xlabel("Year") ax.set_title("PDB entries for EGFR") plt.show()
クエリタンパク質の全PDB IDを取得
次に、pypdb
関数のmake_query
とdo_search
を使って、クエリタンパク質EGFRについての全てのPDB構造を取得します。
search_dict = make_query("P00533") # max_resultsが180以上の場合タイムアウトになる可能性があります。 found_pdb_ids = do_search(search_dict) print("PDB IDs found for query: ") print(found_pdb_ids) print("\nNumber of structures: " + str(len(found_pdb_ids))) # PDB IDs found for query: # ['1IVO', '1M14', '1M17', '1MOX', '1XKK', '1YY9', '1Z9I', '2EB2', '2EB3', '2GS2', '2GS7', '2ITN', '2ITO', '2ITP', '2ITQ', '2ITT', '2ITU', '2ITV', '2ITW', '2ITX', '2ITY', '2ITZ', '2J5E', '2J5F', '2J6M', '2JIT', '2JIU', '2JIV', '2KS1', '2M0B', '2M20', '2N5S', '2RF9', '2RFD', '2RFE', '2RGP', '3B2U', '3B2V', '3BEL', '3BUO', '3C09', '3GOP', '3GT8', '3IKA', '3LZB', '3NJP', '3OB2', '3OP0', '3P0Y', '3PFV', '3POZ', '3QWQ', '3UG1', '3UG2', '3VJN', '3VJO', '3VRP', '3VRR', '3W2O', '3W2P', '3W2Q', '3W2R', '3W2S', '3W32', '3W33', '4G5J', '4G5P', '4HJO', '4I1Z', '4I20', '4I21', '4I22', '4I23', '4I24', '4JQ7', '4JQ8', '4JR3', '4JRV', '4KRL', '4KRM', '4KRO', '4KRP', '4LI5', '4LL0', '4LQM', '4LRM', '4R3P', '4R3R', '4R5S', '4RIW', '4RIX', '4RIY', '4RJ4', '4RJ5', '4RJ6', '4RJ7', '4RJ8', '4TKS', '4UIP', '4UV7', '4WD5', '4WKQ', '4WRG', '4ZAU', '4ZJV', '4ZSE', '5C8K', '5C8M', '5C8N', '5CAL', '5CAN', '5CAO', '5CAP', '5CAQ', '5CAS', '5CAU', '5CAV', '5CNN', '5CNO', '5CZH', '5CZI', '5D41', '5EDP', '5EDQ', '5EDR', '5EM5', '5EM6', '5EM7', '5EM8', '5FED', '5FEE', '5FEQ', '5GMP', '5GNK', '5GTY', '5GTZ', '5HCX', '5HCY', '5HCZ', '5HG5', '5HG7', '5HG8', '5HG9', '5HIB', '5HIC', '5J9Y', '5J9Z', '5JEB', '5LV6', '5SX4', '5SX5', '5U8L', '5UG8', '5UG9', '5UGA', '5UGB', '5UGC', '5UWD', '5WB7', '5WB8', '5X26', '5X27', '5X28', '5X2A', '5X2C', '5X2F', '5X2K', '5XDK', '5XDL', '5XGM', '5XGN', '5XWD', '5Y25', '5Y9T', '5YU9', '5ZTO', '5ZWJ', '6ARU', '6B3S', '6D8E', '6DUK', '6JRJ', '6JRK', '6JRX', '6JWL', '6JX0', '6JX4', '6JXT', '6JZ0', '6P1D', '6P1L', '6P8Q', '6S89', '6S8A', '6S9B', '6S9C', '6S9D', '6V5N', '6V5P', '6V66', '6V6K', '6V6O', '6VH4', '6VHN', '6VHP'] # Number of structures: 205
訳注(05/2020)
RCSB PDBのWebページで「P00533」を検索するとヒットしたエントリーは207でしたが、ここでは205個のエントリーがヒットしています。「2.5.2節:クエリタンパク質のPDBエントリーに関する統計値を取得」においてmax_results=207
とすると途中でとまってしまいましたが、205
では実行可能でした。タイムアウトが原因か、もしくはデータベースの問題かはわかりませんでした。
訳注ここまで
PDBエントリーのメタ情報を取得
describe_pdb
を使って構造のメタ情報を取得します。メタ情報は構造毎に辞書として格納されます。
注:ここではPDB構造のメタ情報を取得するだけで、構造そのもの(3次元座標)はまだ取得しません。
# このデータベースクエリには少し時間がかかるかもしれません。 pdbs = [] for i in found_pdb_ids: pdbs.append(describe_pdb(i)) pdbs[0]
{'relatedPDB': {'@pdbId': '1JL9',
'@details': '1JL9 contains dymeric human EGF molecules.'},
'structureId': '1IVO',
'title': 'Crystal Structure of the Complex of Human Epidermal Growth Factor and Receptor Extracellular Domains.',
'pubmedId': '12297050',
'expMethod': 'X-RAY DIFFRACTION',
'resolution': '3.30',
'keywords': 'TRANSFERASE/SIGNALING PROTEIN',
'nr_entities': '2',
'nr_residues': '1350',
'nr_atoms': '8813',
'deposition_date': '2002-03-28',
'release_date': '2002-10-16',
'last_modification_date': '2011-07-13',
'structure_authors': 'Ogiso, H., Ishitani, R., Nureki, O., Fukai, S., Yamanaka, M., Kim, J.H., Saito, K., Shirouzu, M., Yokoyama, S., RIKEN Structural Genomics/Proteomics Initiative (RSGI)',
'citation_authors': 'Ogiso, H., Ishitani, R., Nureki, O., Fukai, S., Yamanaka, M., Kim, J.H., Saito, K., Inoue, M., Shirouzu, M., Yokoyama, S.',
'status': 'CURRENT'}
PDBエントリーのメタ情報をフィルタリングし並べ替え
取得した情報を関連するPDB構造をフィルタリングして絞り込むために使いたいので、データセットを辞書からより扱いやすいデータフレームに変換します。
pdbs = pd.DataFrame(pdbs) pdbs.head()
relatedPDB | structureId | title | pubmedId | expMethod | resolution | keywords | nr_entities | nr_residues | nr_atoms | deposition_date | release_date | last_modification_date | structure_authors | citation_authors | status | pubmedCentralId | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | {'@pdbId': '1JL9', '@details': '1JL9 contains ... | 1IVO | Crystal Structure of the Complex of Human Epid... | 12297050 | X-RAY DIFFRACTION | 3.30 | TRANSFERASE/SIGNALING PROTEIN | 2 | 1350 | 8813 | 2002-03-28 | 2002-10-16 | 2011-07-13 | Ogiso, H., Ishitani, R., Nureki, O., Fukai, S.... | Ogiso, H., Ishitani, R., Nureki, O., Fukai, S.... | CURRENT | NaN |
1 | {'@pdbId': '1M17', '@details': 'Epidermal Grow... | 1M14 | Tyrosine Kinase Domain from Epidermal Growth F... | 12196540 | X-RAY DIFFRACTION | 2.60 | TRANSFERASE | 1 | 333 | 2452 | 2002-06-17 | 2002-09-04 | 2011-07-13 | Stamos, J., Sliwkowski, M.X., Eigenbrot, C. | Stamos, J., Sliwkowski, M.X., Eigenbrot, C. | CURRENT | NaN |
2 | {'@pdbId': '1M14', '@details': 'Apo-form Epide... | 1M17 | Epidermal Growth Factor Receptor tyrosine kina... | 12196540 | X-RAY DIFFRACTION | 2.60 | TRANSFERASE | 1 | 333 | 2540 | 2002-06-17 | 2002-09-04 | 2011-07-13 | Stamos, J., Sliwkowski, M.X., Eigenbrot, C. | Stamos, J., Sliwkowski, M.X., Eigenbrot, C. | CURRENT | NaN |
3 | [{'@pdbId': '1IGR', '@details': '1IGR contains... | 1MOX | Crystal Structure of Human Epidermal Growth Fa... | 12297049 | X-RAY DIFFRACTION | 2.50 | transferase/growth factor | 2 | 1102 | 8607 | 2002-09-10 | 2003-09-10 | 2011-07-13 | Garrett, T.P.J., McKern, N.M., Lou, M., Ellema... | Garrett, T.P.J., McKern, N.M., Lou, M., Ellema... | CURRENT | NaN |
4 | NaN | 1XKK | EGFR kinase domain complexed with a quinazolin... | 15374980 | X-RAY DIFFRACTION | 2.40 | TRANSFERASE | 1 | 352 | 2299 | 2004-09-29 | 2004-12-07 | 2011-07-13 | Wood, E.R., Truesdale, A.T., McDonald, O.B., Y... | Wood, E.R., Truesdale, A.T., McDonald, O.B., Y... | CURRENT | NaN |
print("Number of PDB structures for EGFR: " + str(len(pdbs))) # Number of PDB structures for EGFR: 205
データセットを以下のクライテリアに沿ってフィルタリングしていきます。
1. 実験手法:X線回折(X-ray diffraction)
X線回折(X-RAY DIFFRACTION
)によって解かれた構造のみを残します。X線回折は最もよく使われている構造決定手法です。
pdbs = pdbs[pdbs.expMethod =="X-RAY DIFFRACTION"] print("Number of PDB structures for EGFR from X-ray: " + str(len(pdbs))) # Number of PDB structures for EGFR from X-ray: 199
2. 構造分解能
分解能が3Å(Angström)以下の構造のみを残します。分解能の値が小さいほど、構造の質が良くなります(=質が良いとは、原子に割り当てられた3次元座標の確からしさが高くなるということです)。3Åよりも下の場合、原子レベルの配置を決めることができ、したがって構造に基づく医薬品デザイン(structure-based drug design)に関連する構造を選択するための閾値として使われることがよくあります。
pdbs_resolution = [float(i) for i in pdbs.resolution.tolist()] pdbs = pdbs[[i <= 3.0 for i in pdbs_resolution]] print("Number of PDB structures for EGFR from X-ray with resolution <= 3.0 Angström: " + str(len(pdbs))) # Number of PDB structures for EGFR from X-ray with resolution <= 3.0 Angström: 164
データセットを構造の分解能で並べ替えます。
pdbs = pdbs.sort_values(["resolution"], ascending=True, na_position='last')
(分解能でソートした)上位のPDB構造を確認します。
pdbs.head()[["structureId", "resolution"]]
structureId | resolution | |
---|---|---|
153 | 5UG9 | 1.33 |
141 | 5HG8 | 1.42 |
152 | 5UG8 | 1.46 |
50 | 3POZ | 1.50 |
56 | 3VRP | 1.52 |
3. リガンドの結合した構造
次のトークトリアルでリガンドベースの組み合わせファーマコフォア(ensemble ligand-based pharmacophore)を作成するので、データフレームから結合するリガンドを含まないPDB構造を取り除きます。pypdb
関数のget_ligands
を使うことでPDB構造に含まれるリガンドの確認と取得ができます。PDBでリガンドと分類されているものにはリガンド、補因子だけでなく溶媒とイオンも含まれています。リガンドの結合した構造だけに絞り込むため、(i) リガンドに分類されている物を含まないすべての構造を取り除き、そして(ii) 分子量(MW)100Da(ダルトン)以上のリガンドを含まない構造をすべて取り除きます。分子量100Daとするのは多くの溶媒とイオンの分子量がそれよりも小さいからです。注:この方法は簡単ですが、全ての溶媒とイオンを除去する方法ではありません。
# データフレームからPDB IDを取得 # pdb_ids = pdbs["structureId"].get_values().tolist() # 訳注: オリジナル資料の上記コードはpandas version 0.25.0以降では使用できなくなっています。 pdb_ids = pdbs["structureId"].to_numpy().tolist()
# 構造の除去 # (i) リガンドを持たないもの # (ii) 分子量(MW)が100 Da(Dalton)以上のリガンドを持たないもの mw_cutoff = 100.0 # 分子量のカットオフ値(単位Da) # このデータベースクエリはしばらく時間がかかるかもしれません removed_pdb_ids = [] for i in pdb_ids: ligand_dict = get_ligands(i) # (i) リガンドがない場合に構造を取り除く if ligand_dict["ligandInfo"] is None: pdb_ids.remove(i) # リストからリガンドのないPDB IDを取り除く removed_pdb_ids.append(i) # リガンドのないPDB IDを格納 # (ii) mw_cutoff以上の分子量のリガンドを一つも持たない構造の除去 else: # リガンド情報の取得 ligs = ligand_dict["ligandInfo"]["ligand"] # 専門的な事項:リガンドが一つの時、辞書型からリスト型へキャストします(続けてリスト内包表記処理するため) if type(ligs) == dict: ligs = [ligs] # リガンドと分類されているもの毎に分子量を取得 mw_list = [float(i["@molecularWeight"]) for i in ligs] # 分子量(MW)100Da(ダルトン)以上のリガンドを含まない構造を除去 if sum([mw > mw_cutoff for mw in mw_list]) == 0: pdb_ids.remove(i) # リストからリガンドのないPDB IDを取り除く removed_pdb_ids.append(i) # リガンドのないPDB IDを格納 print("PDB structures without a ligand (removed from our data set):") print(removed_pdb_ids) # PDB structures without a ligand (removed from our data set): # ['2EB2', '1M14', '2GS2', '2RFE', '4I1Z']
print("Number of structures with ligand: " + str(len(pdb_ids))) # Number of structures with ligand: 159
上位の構造からリガンドのメタ情報を取得
次のトークトリアルでは、最も分解能の高い上位のtop_num
構造からリガンドベースの組み合わせファーマコフォアを作成します。
top_num = 4 # 上位の構造の数 pdb_ids = pdb_ids[0:top_num] pdb_ids # ['5UG9', '5HG8', '5UG8', '3POZ']
get_ligands
を使って上位top_num
個のリガンドについてPDB情報を取得し、csv ファイルとして保存します(リガンド毎の辞書)。
構造が複数のリガンドを含んでいる場合、最も大きなリガンドを選びます。注:これは単純な方法で、タンパク質の結合サイトに結合するリガンドを選択するうえで完璧な方法ではありません。このアプローチではタンパク質に結合する補因子を選んでしまう可能性もあります。ですので、以降のステップで使用する前に、PyMolで自動的に選択した上位のリガンドをチェックしてください。
ligands_list = [] for i in pdb_ids: ligands = get_ligands(i)["ligandInfo"]["ligand"] # 専門的な事項:リガンドが一つの時、辞書型からリスト型へキャストします(続けてリスト内包表記処理するため) if type(ligands) == dict: ligands = [ligands] weight = 0 this_lig = {} # 幾つかリガンドが含まれている場合、最大のものを取得 for lig in ligands: if float(lig["@molecularWeight"]) > weight: this_lig = lig weight = float(lig["@molecularWeight"]) ligands_list.append(this_lig) # フォーマットをデータフレームに変更 ligs = pd.DataFrame(ligands_list) ligs
@structureId | @chemicalID | @type | @molecularWeight | chemicalName | formula | InChI | InChIKey | smiles | |
---|---|---|---|---|---|---|---|---|---|
0 | 5UG9 | 8AM | non-polymer | 445.494 | N-[(3R,4R)-4-fluoro-1-{6-[(3-methoxy-1-methyl-... | C20 H28 F N9 O2 | InChI=1S/C20H28FN9O2/c1-6-15(31)23-13-9-29(7-1... | MJLFLAORJNTNDV-CHWSQXEVSA-N | CCC(=O)N[C@@H]1CN(C[C@H]1F)c2nc(c3c(n2)n(cn3)C... |
1 | 5HG8 | 634 | non-polymer | 377.4 | N-[3-({2-[(1-methyl-1H-pyrazol-4-yl)amino]-7H-... | C19 H19 N7 O2 | InChI=1S/C19H19N7O2/c1-3-16(27)22-12-5-4-6-14(... | YWNHZBNRKJYHTR-UHFFFAOYSA-N | CCC(=O)Nc1cccc(c1)Oc2c3cc[nH]c3nc(n2)Nc4cnn(c4)C |
2 | 5UG8 | 8BP | non-polymer | 415.468 | N-[(3R,4R)-4-fluoro-1-{6-[(1-methyl-1H-pyrazol... | C19 H26 F N9 O | InChI=1S/C19H26FN9O/c1-5-15(30)24-14-9-28(8-13... | CGULPICMFDDQRH-ZIAGYGMSSA-N | CCC(=O)N[C@@H]1CN(C[C@H]1F)c2nc(c3c(n2)n(cn3)C... |
3 | 3POZ | 03P | non-polymer | 547.957 | N-{2-[4-({3-chloro-4-[3-(trifluoromethyl)pheno... | C26 H25 Cl F3 N5 O3 | InChI=1S/C26H25ClF3N5O3/c1-25(2,37)14-22(36)31... | ZYQXEVJIFYIBHZ-UHFFFAOYSA-N | CC(C)(CC(=O)NCCn1ccc2c1c(ncn2)Nc3ccc(c(c3)Cl)O... |
ligs.to_csv('../data/T8/PDB_top_ligands.csv', header=True, index=False, sep='\t')
上位リガンド分子の描画
PandasTools.AddMoleculeColumnToFrame(ligs, 'smiles') Draw.MolsToGridImage(mols=list(ligs.ROMol), legends=list(ligs['@chemicalID']+', '+ligs['@structureId']), molsPerRow=top_num)
タンパク質ーリガンド IDペアの作成
pairs = collections.OrderedDict() for idx, row in ligs.iterrows(): pairs[str(row['@structureId'])] = str(row['@chemicalID']) print(pairs) # OrderedDict([('5UG9', '8AM'), ('5HG8', '634'), ('5UG8', '8BP'), ('3POZ', '03P')])
訳注(2020/05)
OrderedDictは順序つきの辞書です。ここではタンパク質のPDB IDとリガンドのIDをそれぞれDataFrameの@structureId
列と@chemicalID
列から取得し、keyとvalueとしています。
訳注ここまで
PDB構造ファイルの取得
では、pypdb
のget_pdb_file
関数を使ってPDBからPDB構造ファイルを取得します。すなわちタンパク質、リガンド(そして可能な場合には補因子や水分子、そしてイオンといった他の原子あるいは分子)の3次元座標を取得します。
利用可能なファイルフォーマットはpdbとcifで、タンパク質(とリガンド、補因子、水分子、イオン)の原子の3次元座標と、原子間の結合に関する情報とが格納されています。ここではpdbファイルを使用します。
# PDBファイルを取ってきてローカルに保存する for prot, lig in pairs.items(): pdb_file = get_pdb_file(prot, filetype='pdb', compression=False) with open('../data/T8/'+ prot + '.pdb', 'w') as f: f.write(pdb_file)
PDB構造の重ね合わせ
次のトークトリアルでリガンドベースの組み合わせファーマコフォアを作成したいので、全ての構造を互いに3次元で重ね合わせする(アラインメントを取る)必要があります。このタスクのため、分子描画プログラムPyMolを使用します。PyMolはJupyter notebookの中でも使うことができます。PyMolは2つの構造間の原子の距離が最小になるように、一度に2つの構造の重ね合わせをします。
まず最初にコマンドラインからPyMolを立ち上げます(quietモード、すなわちGUIが開かないモード)。
# PyMolをquietモードで立ち上げる pymol.pymol_argv = ['pymol','-qc'] pymol.finish_launching()
重ね合わせのために、リファレンスとなる構造のファイル(動かさないPDB、'target')を選択し、それに対してもう一方の('query')構造ファイルを、PyMolのコマンドcmd.align(query, target)
を使って重ね合わせます。全ての cmd.
コマンドはPyMolのコマンドです。
重ね合わせた構造を新しい座標とともにpdb ファイルとして保存します。また、構造ファイルからリガンドを取り出して別のpdbファイルとして保存し、次のトークトリアルで使用します。
# アラインメントのログをファイルに保存 f = open("../data/T8/alignments.log", "w") # アラインメントの間、動かさない構造と動かす構造とを区別するための変数 immobile_pdb = True refAlignTarget='non' refAlignQuery='non' # 最初のタンパク質に対してタンパク質のアラインメントを取る for prot, lig in pairs.items(): # 動かさない構造(アラインメントのリファレンス構造) if immobile_pdb: target = prot f.write('Immobile target: ' + prot + '\n') # pdbファイルの読み込み(タンパク質とリガンドの複合体) targetFile = cmd.load('../data/T8/' + target + '.pdb') # 精密化したアラインメントのための名前を格納 refAlignTarget='('+target+' within 5 of resn '+lig+')' # 複合体をpdbファイルとして保存 cmd.save('../data/T8/' + target + '_algn.pdb', selection=target) # 選択した名前のリガンドだけを選択 ligObj = cmd.select('ligand', target + ' and resn ' + lig) # 選択したものをpdbファイルとして保存 cmd.save('../data/T8/' + target + '_lig.pdb', selection='ligand', format='pdb') # リガンドの選択を削除 cmd.delete(ligObj) # targetが選択されました。 immobile_pdb = False # 動かす構造(リファレンス構造に対して重ね合わせる) else: query = prot f.write('-- align %s to %s \n' %(query, target)) # pdbファイルの読み込み(タンパク質のリガンドの複合体) queryFile = cmd.load('../data/T8/' + query + '.pdb') # 結合サイトに焦点をあてて構造(タンパク質)を重ね合わせる refAlignQuery= '('+query+' within 5 of resn '+lig+')' values = cmd.align(refAlignQuery, refAlignTarget) # 構造の重ね合わせができない場合(即ち、RMSD > 5A の場合)、スキップする if values[0] > 5: f.write('--- bad alignment: skip structure\n') else: # 複合体をpdbファイルとして保存 cmd.save('../data/T8/' + query + '_algn.pdb', selection=query) # リガンドのみを選択 ligObj = cmd.select('ligand', query + ' and resn ' + lig) # 選択したものをpdbファイルとして保存 cmd.save('../data/T8/' + query + '_lig.pdb', selection='ligand', format='pdb') # リガンドの選択を削除 cmd.delete(ligObj) # "query"の選択を削除 cmd.delete(query) # "target"の選択を削除 cmd.delete(target) f.close()
訳注(2020/05)
Pymol Wikiを参照して上のセルの補足をします。
アラインメントのコマンドcmd.align
に渡す引数のrefAlignTarget
とrefAlignQuery
はどちらも( PDB_ID within 5 resn of Ligand_ID)
という文字列になる変数です。これは「タンパク質(PDB_ID)のうち、残基名(resn, residune name)で指定されたリガンド(Ligand_ID)の5Å以内に入る部分」を指し示すようです。リガンドの周囲を指定してアラインメントをとっているので「結合サイトに焦点をあてて重ね合わせる」ことになるようです。
また、cmd.align
の戻り値は7つあり、一番最初(index 0)が「精密化した後のRMSD」となっているため、if values[0] > 5
という条件文になっているようです(PyMol Wiki: Align)。
訳注ここまで
PyMolアプリケーションを終了します。
# PyMolを終了
pymol.cmd.quit()
全てのリガンドpdbファイルが存在していることを確認します。もしファイルが欠けていたなら、手作業でPyMOLでタンパク質ーリガンド複合体構造をチェックしてください。
mol_files = [] for file in glob.glob("../data/T8/*_lig.pdb"): mol_files.append(file) mol_files # ['../data/T8/5UG8_lig.pdb', # '../data/T8/3POZ_lig.pdb', # '../data/T8/5UG9_lig.pdb', # '../data/T8/5HG8_lig.pdb']
ディスカッション
このトークトリアルでは、PDBからタンパク質とリガンドのメタ情報と構造情報を取得する方法を学びました。X線構造だけを残し、分解能とリガンドがあるかどうかによって絞りこみました。最終的な目標は、重ね合わせたリガンドのセットを使って、次のトークトリアルでリガンドベースの組み合わせファーマコフォアを作成することです。
ファーマコフォアモデリングのためにリガンドについての情報量を多くするため、PDB構造を分解能でフィルタリングするだけでなく、リガンドの多様性をチェックすること(類似性による分子クラスタリングについてはトークトリアル 5を参照)と、リガンドの生理活性をチェックする(即ち、活性があるリガンドだけを含める)ことをお勧めします。
クイズ
- Protein Data Bankの含むデータの種類を要約してください。
- 構造の分解能が何を表し、このトークトリアルではどうやって、また、なぜ分解能についてフィルタリングしたのかを説明してください。
- 構造の重ね合わせ(アラインメントを取ること)とはどういう意味か説明し、このトークとリアルで実施したアラインメントについて議論してください。
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回、 取得したPDBの構造を絞り込むうえで分解能が用いられていましたが、登録されている構造の正確さについては色々と注意が必要だそうです。
@torusengoku先生の記事: X線結晶構造中の低分子や水の妥当性評価 で驚きの事例を紹介してくださっています。
PDBからの複数の構造について、KNIMEを活用する方法を 或る化みす途のブログ さんが紹介してくださっています。
* Dear Pymol Beginner x KNIME -KNIMEを使ってPDBファイルをまとめてDLしてみたり-
PyMol 機能ありすぎてよくわからないのですが、とりあえず pymol-book が素晴らしいので読んでみようと思います(知らなかった・・・)。
TeachOpenCADD トピック7 〜リガンドベーススクリーニング〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル7をもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!
以下、訳。
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
トークトリアル 7
リガンドベーススクリーニング:機械学習
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Jan Philipp Albrecht and Jacob Gora
このトークトリアルの目的
利用可能なデータソースがより大きくなるにつれ、医薬品探索、特にリガンドベースのバーチャルスクリーニングにおいて、機械学習(machine learning、ML)が盛り上がってきました。このトークトリアルでは、私たちの標的とするターゲット分子(EGFR)に対して新規な化合物の活性を予測するために、様々な教師あり機械学習(supervised ML)アルゴリズムを使用する方法について学習します。
学習の目標
理論
実践
- EGFR阻害剤となる可能性のある化合物を見つけるための、機械学習に基づくスクリーニングパイプラインの構築と評価
レファレンス
- RdKitフィンガープリント:
例)2012年RDKit UGMでのG. Landrumによるプレゼンテーションを参照してください - 機械学習:
- ランダムフォレスト(Random forest、RF): http://ect.bell-labs.com/who/tkh/publications/papers/odt.pdf
- サポートベクトルマシン(Support vector machines、SVM): https://link.springer.com/article/10.1007%2FBF00994018
- 人工ニューラルネットワーク(Artificial neural networks、ANN): https://www.frontiersin.org/research-topics/4817/artificial-neural-networks-as-models-of-neural-information-processing
- パフォーマンス:
- J. Med. Chem., 2017, 60, 474−485 に基づくB. MergetによるGitHubのNotebook も参照してください。
理論
機械学習をうまく適用するために必要となるものは、巨大な化合物データセット、化合物のエンコーディング、データセットの化合物ごとのラベル、そしてモデルを訓練するための機械学習アルゴリズムです。これらによって、新しい化合物の予測を行うことができます。
データの準備
機械学習のために、化合物を特徴のリストに変換する必要があります。化合物の表現方法として、フィンガープリントがよく使われます。
このトークトリアルで使うフィンガープリントはRDKitに実装されているものです(より詳細な情報はG. Landrumのプレゼンテーション で見ることができます)。
- maccs: MACCS キー 166ビットの構造キー記述子(structural key descriptors)で、各ビットはSMARTSパターンに関連づけられています。
- ecfp4とecfp6: Extended-Connectivity Fingerprints (ECFPs) は化合物の特徴づけ、類似性検索、そして構造活性モデリングのためにデザインされた円形トポロジカルフィンガープリント(circular topological fingerprints)です。ECFPの最も重要なパラメータは最大直径とフィンガープリントの長さです。直径(diameter)と呼ばれるパラメータは、各原子について考慮する円形状の近傍の最大の直径を明示します。ここでは2つの直径、4と6があります。長さ(length)のパラメータはビット文字列表現の長さを明示します。デフォルトの長さは2048です。
- torsion: Torsion Fingerprint Deviation (TFD) は、クエリ分子と、非環状の結合(acyclic bond)と環系(ring system)を考慮して生成した配座群からトーションフィンガープリント(Torsion Fingerprint)を抽出、重み付けし、そして比較します。
- rdk5: rdk5は経路に基づくフィンガープリント(path based fingerprint)です。パスフィンガープリントは、分子グラフの全ての直線的なフラグメントを与えられたサイズまで網羅的に列挙し、これらのフラグメントを固定長ビットベクトルにハッシュ化することで生成されます。
機械学習 (Machine Learning、ML)
機械学習は次の用途に使うことができます(scikit-learn pageも参照してください)。
- (教師あり)分類:あるオブジェクト(object)がどのカテゴリーに属しているかを明らかにする(最近傍、ナイーブベイズ、ランダムフォレスト、サポートベクトルマシン、・・・)
- 回帰:オブジェクトに関連づけられた連続値の特徴量の予測
- (教師なし)クラスタリング:類似のオブジェクトを各セットに自動的にグループ化( トークトリアル5 参照)
教師あり学習
学習アルゴリズムは訓練データの中のパターンを見つけることでルールを作成します。
- ランダムフォレスト(Random Forest、RF):多数の決定木(decision tree)により平均化された予測を出します。
- サポートベクトルマシン(Support Vector Machines (SVM)): SVMはカーネルトリックと呼ばれる入力を高次元の特徴空間に暗にマッピングする手法を使って、非線形の分類を効率的に実行することができます。オブジェクトの関数としてマージンを最大化するというアイデアに基づいた分類です。
- 人工ニューラルネットワーク(Artificial neural networks、ANNs):ANNは生物学的な脳の神経をゆるくモデル化した人工神経(artificial neuron)と呼ばれる連結されたユニットあるいはノードの集まりに基づいています。生物学的な脳のシナプスの様に、各結合は一つの人工神経からもう一方へとシグナルを伝達することができます。シグナルを受け取った人工神経は、そのシグナルを処理し、結合する別の人工神経へとシグナルを送ります。(図はWikipediaより)
妥当性検査手法(Validation strategy):k分割交差検証
このモデル妥当性評価テクニックは反復的な方法でデータセットを2つのグループに分けます。
この検証の目標は過学習(over-fitting)として知られる問題に注意するために、モデルが一度も見たことのないデータを予測する性能をテストし、モデルの汎化性能を評価することです。
性能評価
- 感度(Sensitivity)、または真陽性(true positive rate): TPR = TP/(FN+TP)
- 特異度(Specificity)、または真陰性(true negative rate): TNR = TN/(FP + TN)
- 精度(Accuracy)、または真度(trueness): ACC = (TP + TN)/(TP + TN + FP + FN)
- ROC曲線、受信者操作特性曲線(receiver operating characteristic curve、ROC-curve):
- 分類器(classifier)の診断性能を図解するプロット図
- 感度を特異度に対してプロットします。
- AUC、ROC曲線の下の面積(Area Under the roc Curve):
- 分類器がランダムに選択した正例を負例よりも高くランクづけする確率を表現します。
- 0から1の値をとり、高い方がより良い値です。
# インポート文 # 一般的な項目: import pandas as pd import numpy as np # RDKit: from rdkit import Chem from rdkit.Chem import RDKFingerprint from rdkit.Chem.AllChem import GetMorganFingerprintAsBitVect from rdkit.Chem.AllChem import GetHashedTopologicalTorsionFingerprintAsBitVect from rdkit.Chem import MACCSkeys from rdkit.DataStructs import ConvertToNumpyArray # sklearn: from sklearn.ensemble import RandomForestClassifier from sklearn import svm from sklearn.neural_network import MLPClassifier from sklearn.model_selection import KFold from sklearn.metrics import auc from sklearn.metrics import accuracy_score from sklearn.metrics import recall_score from sklearn.metrics import roc_curve # from sklearn.manifold import MDS # matplotlib: import matplotlib.ticker as ticker import matplotlib.pyplot as plt import matplotlib.cm as cm # seaborn: import seaborn as sns
データの準備
ここではEGFR(Epidermal growth factor receptor)キナーゼのデータを扱います。
データを取り扱う前にまず、作業を行うデータフレームの作成に役立つ2つの関数を定義します。
最初の関数は分子のフィンガープリントを計算するもので calculate_fp
と名付けます。フィンガープリントの選択肢は次です。
- maccs
- ecfp4 と ecfp6
- torsion
- rdk5
def calculate_fp(mol, method='maccs', n_bits=2048): # mol = Chem molecule object # 与えられたビットの数と手法で化合物のフィンガープリントを計算する関数 if method == 'maccs': return MACCSkeys.GenMACCSKeys(mol) if method == 'ecfp4': return GetMorganFingerprintAsBitVect(mol, 2, nBits=n_bits, useFeatures=False) if method == 'ecfp6': return GetMorganFingerprintAsBitVect(mol, 3, nBits=n_bits, useFeatures=False) if method == 'torsion': return GetHashedTopologicalTorsionFingerprintAsBitVect(mol, nBits=n_bits) if method == 'rdk5': return RDKFingerprint(mol, maxPath=5, fpSize=1024, nBitsPerHash=2)
2番目の関数を使うことで次の列が追加されたデータフレームを作成することができます:
- moleculeオブジェクトとしての分子(SMILES文字列から作成)
- 選択した手法のPythonオブジェクトとしてのフィンガープリント(ここではMACCs)
- 選択したフィンガープリント手法のビットベクトルのバイナリ表記
したがって、パラメータが2つあります。一つは、正しいSMILES表記を格納する「smiles」という名前の列と、フィンガープリントの長さを保持するデータフレームです。2番目の引数はMACCSフィンガープリント以外の手法を選択した場合に使用します。
def create_mol(df_l, n_bits): # SMILES文字列から分子を構築 # mol列を生成:molオブジェクトを返し、失敗した場合はNoneを返す df_l['mol'] = df_l.smiles.apply(Chem.MolFromSmiles) # 化合物フィンガープリントをフィンガープリントオブジェクトとして格納する列を生成 df_l['bv'] = df_l.mol.apply( # 各化合物についてラムダ関数「calculate_fp」を適用 lambda x: calculate_fp(x, 'maccs', n_bits) ) # np.array をフィンガープリントのビットベクトルを保持するために割り当て(np = numpy) df_l['np_bv'] = np.zeros((len(df_l), df_l['bv'][0].GetNumBits())).tolist() df_l.np_bv = df_l.np_bv.apply(np.array) # フィンガープリントオブジェクトをNumpuArrayに変換しnp_bvに格納 df_l.apply(lambda x: ConvertToNumpyArray(x.bv, x.np_bv), axis=1)
データの読み込み
さて、データを読み込み、実際の作業を始めましょう。トークトリアル 2で取得したcsvファイルを、重要な列とともにデータフレームに読み込みます:
- CHEMBL-ID
- 対応する化合物のSMILES値
- pIC50
# 先のトークトリアルのデータを読み込む df = pd.read_csv('../data/T2/EGFR_compounds_lipinski.csv', delimiter=';', index_col=0) # 先頭行をみる print(df.shape) print(df.info()) # (4523, 10) # <class 'pandas.core.frame.DataFrame'> # Int64Index: 4523 entries, 0 to 5427 # Data columns (total 10 columns): # molecule_chembl_id 4523 non-null object # units 4523 non-null object # IC50 4523 non-null float64 # smiles 4523 non-null object # pIC50 4523 non-null float64 # MW 4523 non-null float64 # HBA 4523 non-null int64 # HBD 4523 non-null int64 # LogP 4523 non-null float64 # rule_of_five_conform 4523 non-null object # dtypes: float64(4), int64(2), object(4) # memory usage: 388.7+ KB # None
データを分類
各化合物を活性あるいは不活性として分類する必要があるので、pIC50値を使います。 * pIC50 = -log10(IC50) * IC50はin-vitroの阻害が50%になるモル濃度(mol/L)を表します。 * pIC50のデータを離散化(活性と不活性に分類)するのによく用いられるカットオフ値は6.3で、私たちの実験でもこの値を使います。 * 文献上ではpIC50値5から7までの範囲でいくつか他にも活性のカットオフが提案されており、データポイントをとらない除外領域を定義しているものさえある、ということを言及しておきます。
これで、上で定義した関数を使って、化合物とそのフィンガープリントを生成し、どの化合物が活性で、どれが活性でないかを明示することができます。
# 不必要な列の削除 df_new=df.drop(['units', 'IC50'], axis=1) # SMILESから化合物とフィンガープリントを作成 create_mol(df_new, 2048) # 活性を格納する列を追加 df_new['active'] = np.zeros(len(df_new)) # pIC50が6.3より大きい全ての化合物に印をつける df_new.loc[df_new[df_new.pIC50 >= 6.3].index, 'active'] = 1.0 print('actives: %d, inactives: %d' % (df_new.active.sum(), len(df_new)-df_new.active.sum())) # actives: 2555, inactives: 1968
df_new.head(3)
molecule_chembl_id | smiles | pIC50 | MW | HBA | HBD | LogP | rule_of_five_conform | mol | bv | np_bv | active | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | CHEMBL63786 | Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 | 11.522879 | 349.021459 | 3 | 1 | 5.2891 | yes | <rdkit.Chem.rdchem.Mol object at 0x12a734a30> | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ... | [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ... | 1.0 |
1 | CHEMBL53711 | CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 | 11.221849 | 343.043258 | 5 | 1 | 3.5969 | yes | <rdkit.Chem.rdchem.Mol object at 0x12a7359e0> | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ... | [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ... | 1.0 |
2 | CHEMBL35820 | CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC | 11.221849 | 387.058239 | 5 | 1 | 4.9333 | yes | <rdkit.Chem.rdchem.Mol object at 0x12a735710> | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ... | [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ... | 1.0 |
機械学習
以下では、我々の化合物を分類するためにいくつかの機械学習のアプローチを試します。テストするアプローチは以下です。
* ランダムフォレスト(Random Forest、RF)
* サポートベクトルマシン(Support Vector Machines、SVM)
* 人工ニューラルネットワーク(Artificial Neural Networks、ANNs)
加えて、結果について評価し議論します。始める前に、 crossvalidation
という名前の関数を定義します。この関数は交差検証(クロスバリデーション)の手順を実行し、精度と感度、特異度といった指標を返します。
目標は過学習(over-fitting)として知られる問題に注意するために、モデルが一度も見たことのないデータを予測する性能をテストし、モデルの汎化性能を評価することです。
# 交差検証ループのための関数 def crossvalidation(model_l, df_l, n_folds=10): # 選択したモデル、データフレームそして分割(fold)数を与えると、この関数は交差検証を実行し、 # 予測の精度、感度、特異度を返し、また、各分割(fold)について偽陽性率(fpr)、真陽性率(tpr)、ROC曲線も返します。 # resultsベクトルを空にする results = [] # k-分割交差検証のためにインデックスをシャッフルする kf = KFold(n_splits=n_folds, shuffle=True) # 各データポイントについてラベルを-1で初期化 labels = -1 * np.ones(len(df_l)) # 分割(fold)数についてループを回す for train_index, test_index in kf.split(df_l): # トレーニング # ビットベクトルとラベルをリストに変換 train_x = df_l.iloc[train_index].bv.tolist() train_y = df_l.iloc[train_index].active.tolist() # モデルを適合させる model_l.fit(train_x, train_y) # テスト # ビットベクトルとラベルをリストに変換 test_x = df_l.iloc[test_index].bv.tolist() test_y = df_l.iloc[test_index].active.tolist() # テストセットについて予測 prediction_prob = model_l.predict_proba(test_x)[:, 1] # 各分割(fold)について予測したラベルを保存 labels[test_index] = model_l.predict(test_x) # 性能評価(パフォーマンス) # 各分割(fold)についての偽陽性率(fpr)、真陽性率(tpr)、ROC曲線 fpr_l, tpr_l, _ = roc_curve(test_y, prediction_prob) roc_auc_l = auc(fpr_l, tpr_l) # 結果に追加 results.append((fpr_l, tpr_l, roc_auc_l)) # 全体の精度、感度、特異度を取得 y = df_l.active.tolist() acc = accuracy_score(df_l.active.tolist(), labels) sens = recall_score(df_l.active.tolist(), labels) spec = (acc * len(y) - sens * sum(y)) / (len(y) - sum(y)) return acc, sens, spec, results
もちろん、私たちはモデルの質を評価したいと思っています。したがって予測の精度と感度、特異度を知りたいと思います。加えていわゆるROC曲線に焦点を当てます。
コードを明確にし理解しやすくするために、結果をプロットする簡単な関数を作ります。
以下の点にひとまず焦点をあてます:
- 感度
- 特異度
- 精度
- ROC曲線とAUC
np.linspace(0.1, 1.0, 10) # array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
def print_results(acc, sens, spec, stat_res, main_text, file_name, plot_figure=1): plt.figure(plot_figure, figsize=(7, 7)) cmap = cm.get_cmap('Blues') colors = [cmap(i) for i in np.linspace(0.3, 1.0, 10)] #colors = ["#3465A4"] for i, (fpr_l, tpr_l, roc_auc_l) in enumerate(stat_res): plt.plot(fpr_l, tpr_l, label='AUC CV$_{0}$ = {1:0.2f}'.format(str(i),roc_auc_l), lw=2, color=colors[i]) plt.xlim([-0.05, 1.05]) plt.ylim([-0.05, 1.05]) plt.plot([0, 1], [0, 1], linestyle='--', label='Random', lw=2, color="black") # Random curve plt.xlabel('False positive rate', size=24) plt.ylabel('True positive rate', size=24) plt.title(main_text, size=24) plt.tick_params(labelsize=16) plt.legend(fontsize=16) # Sプロットを保存 - テキストボックスを含めるためにbbox_inchesを使います。 # https://stackoverflow.com/questions/44642082/text-or-legend-cut-from-matplotlib-figure-on-savefig?rq=1 plt.savefig("../data/T7/" + file_name, dpi=300, bbox_inches="tight", transparent=True) plt.show() # AUCの平均を計算し出力します m_auc = np.mean([elem[2] for elem in r[3]]) print('Mean AUC: {}'.format(m_auc)) # 全体の精度、感度、特異度を表示します print('Sensitivity: {}\nAccuracy: {}\nSpecificity: {}\n'.format(acc, sens, spec)) print('\n')
ランダムフォレスト分類器
さて、ランダムフォレスト分類器からはじめます。まず最初にパラメータを設定します。その後で、モデルのクロスバリデーションを実行し、結果をプロットします。
# ランダムフォレストのパラメーターを設定 param = {'max_features': 'auto', 'n_estimators': 2000, 'criterion': 'entropy', 'min_samples_leaf': 1} modelRf = RandomForestClassifier(**param) # 10分割で交差検証を実行 r = crossvalidation(modelRf, df_new, 10)
# AUCの結果をプロット # r は精度(acc)、感度(sens)、特異度(spec)と結果を保持しています print_results(r[0], r[1], r[2], r[3], 'Random forest ROC curves', 'rf_roc.png', 3) # Mean AUC: 0.8997931471381779 # Sensitivity: 0.83440194561132 # Accuracy: 0.8806262230919765 # Specificity: 0.7743902439024393
我々のモデルは全ての指標について非常に良い値を示し、予測性能があるように見えます。
サポートベクトル分類器
ここでは動径基底関数カーネル(Radial-basis function kernel(あるいは2乗指数カーネル(squared-exponentioal kernel)とも呼ばれます)によるサポートベクトルマシンを学習させます。より詳しい情報は sklearn RBF kernel を参照してください。
# モデルを明示 modelSvm = svm.SVC(kernel='rbf', C=1, gamma=0.1, probability=True) # 10分割交差検証を実行 r = crossvalidation(modelSvm, df_new, 10)
# 結果をプロット print_results(r[0], r[1], r[2], r[3], 'SVM$(rbf kernel)$ $C=1$ $\gamma=0.1$ ROC curves', 'svm_roc.png', 3) # Mean AUC: 0.8969024290217835 # Sensitivity: 0.8441300022109219 # Accuracy: 0.9056751467710372 # Specificity: 0.7642276422764228
ニューラルネットワーク分類器
ここで試す最後のアプローチはニューラルネットワークモデルです。各層5つのニューロンをもつ3層のMLPClassifier (多層パーセプトロン分類器、Multi-layer Perceptron classifier)を学習させます。早期終了(early stopping)が明示的にFALSEに設定されていることに気づかれたかもしれません。これまでと同様、交差検証を実行し結果をプロットします。MLPに関するより詳細な情報はsklearn MLPClassifierを参照してください。
# モデルを明示、デフォルトの活性化関数:relu modelClf = MLPClassifier(solver='adam', alpha=1e-5, hidden_layer_sizes=(5, 3), random_state=1, early_stopping=False) # 10分割交差検証を実行 r = crossvalidation(modelClf, df_new, 10)
# 結果をプロット print_results(r[0], r[1], r[2], r[3], 'MLPClassifier ROC curves', 'mlp_roc.png', 3) # Mean AUC: 0.8686086737582928 # Sensitivity: 0.8032279460535043 # Accuracy: 0.8583170254403131 # Specificity: 0.7317073170731707
ディスカッション
我々のデータッセットで最もパフォーマンスがよかったのはどのモデルで、それは何故でしょうか?
- データセットに対して、3つのモデルはすべて(非常に)良いパフォーマンスでした。最も良いモデルはランダムフォレストとサポートベクトルマシンで、AUCの平均は90%でした。ニューラルネットワークはすこし低い結果で、AUCの平均は87%でした。(スクリプトを再度実行すると少し異なる結果となる可能性があることに注意してください。)
- ランダムフォレストとサポートベクトルマシンのモデルが最も良いパフォーマンスを示したのにはいくつか理由があるかもしれません。我々のデータセットは、単純な決定木用の分類器(tree-like decisions)、あるいは動径基底関数でそれぞれ簡単に活性/不活性を分類できるものだったのかもしれません。したがって、この分類を行ううえでフィンガープリントにあまり複雑なパターンがなかったのかもしれません。
- 人工ニューラルネットワーク(ANN)のパフォーマンスが少し悪かった理由としては、モデルを訓練するためのデータが単純に少なすぎたのかもしれません。
- 加えて、モデルの評価のためには他の外部検証用データセット(validation set)を持っておくことを常におすすめします。
MACCSを用いることは正しい選択だったのでしょうか?
- 明らかに、MACCSは分類が可能かどうかみるために、モデルの訓練と評価を行う出発点として良いものです。
- ですが、MACCSフィンガープリント(166ビット)は他のフィンガープリント(2048ビット)と比較してかなり短いので、異なるフィンガープリントを試し、検証のプロセスを繰り返した方が良いです。
ここから次にどちらに進めば良いか?
- いくつかのモデルをうまく学習させることができました。
- 次のステップは、新奇なEGFR阻害剤となりうる可能性がある化合物を予測するために、これらのモデルを使って未知のスクリーニングデータセットの分類を行うことです。
- 巨大なスクリーニングデータセットの例として、例えば700万以上の化合物からなるMolPortがあります。
- 我々のモデルを使ってMolPortの化合物をランク付けし、活性がある可能性が最も高いと予測された化合物群についてさらに研究を進めることができます。
- そのような適用例としては、S. RinikerとG. Landrumによって開発された TDT Tutorial も参照してみてください。新しい抗マラリア薬を見つけるために、融合モデルを訓練し eMolecules のスクリーニングを行なっています。
クイズ
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。
いつも引用させていただいている「化学の新しいカタチ」さんの記事より
* 化学の新しいカタチ:RDKitとランダムフォレスト
* 化学の新しいカタチ:交差検証を用いてElastic Netを化合物の溶解度データに対して最適化
化学分野と機械学習という点では明治大学金子先生の研究室HPが非常に充実していると思います。 取り上げられた内容をまとめてくださっているページもあるので、そちらを参照いただいて関連する記事にアクセスするのが一番早いかと思います。
ざっと見たところ今回の内容と関連するところとしては
- ランダムフォレスト(Random Forest, RF)~アンサンブル学習で決定木の推定性能を向上!~
- サポートベクターマシン(Support Vector Machine, SVM)~優秀な(非線形)判別関数~
- 回帰モデル・クラス分類モデルを評価・比較するためのモデルの検証 (Model validation)
などがありそうです。素敵!
勝手にいろいろと引用させていただいておりますので問題があれば削除します。
TeachOpenCADD トピック6 〜最大共通部分構造〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル6をもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下、訳。
トークトリアル 6
最大共通部分構造(Maximum common substructure)
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Oliver Nagel
このトークトリアルの目的
サイズの大きい化学データのクラスタリングとクラス分類は、医薬品探索の化学を利用する広範な領域において、研究を導き、分析し、知識を発見することにとって不可欠です。
一つ前のトークトリアルで、化合物をグループ化する方法(クラスタリング)を学び、一つのクラスターに含まれる化合物がお互いに似通っており、共通の骨格(scaffold)を共有していることを見つけました。視覚的に調べることに加えて、ここでは化合物セットが共通してもつ最大の部分構造を計算する方法を学びます。
学習の目標
- 最大共通部分構造とは何か(Maximum Common Substructure、MCS)?
- 化合物のグループのMCSを計算するにはどうすれば良いか?
理論
- 化合物セットの最大共通部分構造をみつけることについての導入
- FMCSアルゴリズムの詳細な説明
実践
レファレンス
- Dalke A, Hastings J., FMCS: a novel algorithm for the multiple MCS problem. J. Cheminf. 2013; 5(Suppl 1): O6
- Raymond JW., Willet P., Maximum common subgraph isomorphism algorithms for the matching of chemical structures. J Comput Aided Mol Des. 2002 Jul;16(7):521-33
- アルゴリズムに関する情報の書かれたDalkeのWebサイト: http://dalkescientific.com/writings/diary/archive/2012/05/12/mcs_background.html
- MCSに関するRDKitクックブックのドキュメンテーション: http://www.rdkit.org/docs/Cookbook.html#using-custom-mcs-atom-types
理論
導入
最大共通部分構造(maximum common structure、MCS)は2つあるいはそれ以上の対象化合物に含まれる最大の部分構造として定義されています。
* MCSをみつけること = グラフ内の最大共通パターンマッチング問題(maximum common subgraph isomorphism problem)
* ケモインフォマティクスの分野で多くの適用用途があります。例えば、類似性検索(similarity search)、階層クラスタリング(hierarchical clustering)、分子アラインメント(molecule alignment)など。
* 利点:
* 直観的 → 対象化合物に共有されている構造は重要である可能性が高い
* ありうる活性のパターンへの洞察を与える
* 単純に部分構造をハイライトすることで簡単に可視化できる
MCSアルゴリズムの詳細 (レビューを参照してください: J Comput Aided Mol Des. 2002 Jul;16(7):521-33)
* 2つあるいはそれ以上のグラフ間のMCSを決めることはNP-完全問題です。
* 厳密に決定するためのアルゴリズムと近似アルゴリズムのどちらもあります。
* 厳密:最大クリーク(Maximum-clique)、後戻り法(バックトラック法、backtracking)、動的計画法(dynamic programming)
* 近似:遺伝的アルゴリズム(Genetic algorithm)、組合せ最適化(combinatorial optimization)、フラグメントストレージ(fragment storage)、・・・
* 問題分割:分子グラフの単純化
FMCS アルゴリズム:
* MCS問題をグラフパターンマッチング問題(graph isomorphism problem)としてモデル化します。
* サブグラフの数え上げとサブグラフのパターンマッチングテスト(subgraph isomorphism testing)に基づきます。
FMCSアルゴリズムの詳細な説明
J. Cheminf. 2013; 5(Suppl 1): O6 とrdkit fmcs documentaion の各項目で説明されています。
単純化したアルゴリズムの説明
best_substructure = None pick one structure in the set as query, all other as targets # 化合物セットから一つをクエリとして選び、残り全てをターゲットにする for each substructure in the query: # クエリの各部分構造について convert into a SMARTS string based on the desired match properties # マッチングを希望する特徴に基づきSMARTS文字列に変換する if SMARTS pattern exists in all of the targets: # もしターゲットのすべてにそのSMARTSパターンが存在していれば then it is a common substructure # それが共通部分構造となります keep track of the maximum of such substructure # そのような部分構造の最大となるものを順に辿っていきます。
通常、この単純なアプローチだけではとても長い時間がかかりますが、実行時間を加速するためのトリックがいくつかあります。
結合除去(Bond elimination)
- MCSの部分とはなり得ない結合の除去
- 各入力構造に原子タイプと結合タイプの情報が存在している必要があります
- 結合タイプ:最初の原子、結合、2つめの原子のSMATSからなる文字列
- 入力構造全てに存在していない全ての結合タイプを排除し、それぞれの辺(edge)を削除
- 結果:全ての原子の情報をもちつつ、辺(edge(結合))の数がより少ないフラグメント化された構造
最大フラグメントのうち最小の構造をクエリとして用いる
- ヒューリスティックなアプローチ:
- 各入力構造の最大のフラグメントを見つける
- 最大フラグメントの結合の数によって入力構造を昇順に並べ換える
- 原子の数との紐付けを解き、入力の順を代わりのものとして使う
- 最大フラグメントのうち最小の構造をクエリ構造とする
- 他の入力構造由来の構造がターゲットとなります
フラグメントのサブグラフを列挙するために、幅優先探索(breadth-first search (BFS))と優先度付きキュー(priority que)を使う
- いわゆる種(シード、seed)を成長させることに基づく列挙(enumeration)
- シード:現在のサブグラフの原子/結合と、除外するセット(成長には使えない結合)
- 冗長性を防ぐために:
- 最初のシードはフラグメントの最初の結合で、フラグメント全体のサイズにまで成長する可能性があります
- 2番目のシードは2番目の結合で、最初の結合を使ったものからは排除されます
- 3番目のシードは3番目の結合から始まり、最初と2番目を使用したものからは排除されます
- ・・・
- シードはシードにつながっている結合に沿って成長します(除外セットや既にシードに含まれているものは含みません)
- 各ステップですべての成長の可能性が考慮されます
- 例、拡張の方向としてN個の可能な結合がある時、2のN-1乗個の新しいシードがキューにつけ加えられる
- サブグラフに加える新しい結合がなくなった時点で列挙は終わります(シードがキューから除外されます)
- 最大のシードが最初に処理されます
他の全ての構造に含まれていないシードを取り除く
- 各成長の段階で → 新しいシードが他の全ての構造に存在しているかどうかチェックする
- 存在していない場合:キューからシードを取り除く
十分な成長のポテンシャルがないシードを取り除く
- 除外リストと拡張の可能性がある辺(エッジ)から成長ポテンシャルを評価します
- 現在の最も良いサブグラフよりもポテンシャルが小さい場合 -> シードをキューから取り除きます
これらのアプローチを使うことで、最大共通部分構造に相当する最大のサブグラフを順に辿っていくことことは普通の問題となります。
実践
# 必要なモジュールをインポート import matplotlib.pyplot as plt from rdkit import Chem from rdkit.Chem import AllChem from collections import defaultdict from rdkit.Chem import Draw from rdkit.Chem.Draw import IPythonConsole from rdkit.Chem import rdFMCS import pandas as pd import numpy from rdkit.Chem import PandasTools IPythonConsole.ipython_useSVG=True
化合物を読み込み描画
suppl = Chem.SDMolSupplier('../data/T5/molSet_largestCluster.sdf') mols = [x for x in suppl] print('Set with %d molecules loaded'%(len(mols))) # 最初の12個の化合物を表示 #Draw.MolsToGridImage(mols[:12], legends=[x.GetProp("_Name") for x in mols[:12]], molsPerRow=4) Draw.MolsToGridImage(mols[:12], molsPerRow=4) # Set with 143 molecules loaded
様々なパラメーター入力値でRDKitに実装されているFMCSアルゴリズムを実行
デフォルト値
最も単純な場合、化合物のリストだけをパラメーターとして与えます。
mcs1 = rdFMCS.FindMCS(mols) print('SMARTS string: %s'%mcs1.smartsString) # 構造の描画 m1 = Chem.MolFromSmarts(mcs1.smartsString) Draw.MolToImage(m1) # SMARTS string: [#6]-[#6](=[#8])-[#7]-[#6]1:[#6]:[#6]2:[#6](-[#7]-[#6]3:[#6]:[#6]:[#6]:[#6]:[#6]:3):[#7]:[#6]:[#7]:[#6]:2:[#6]:[#6]:1-[#8]-[#6]-[#6]
MCSをハイライトして化合物を描画するためのヘルパー関数を定義します。
# クエリー化合物群のMCSをハイライトします def highlightMolecules(cur_mols, cur_mcs, num, label=True): pattern = Chem.MolFromSmarts(cur_mcs.smartsString) matching = [cur_mols[i].GetSubstructMatch(pattern) for i in range(0,len(cur_mols))] if label==True: return Draw.MolsToGridImage(cur_mols[:num], legends=[x.GetProp("_Name") for x in mols[:num]], molsPerRow=3, highlightAtomLists = matching[:num], subImgSize=(300,250)) else: return Draw.MolsToGridImage(cur_mols[:num], molsPerRow=3, highlightAtomLists = matching[:num], subImgSize=(300,250))
highlightMolecules(mols, mcs1, 12)
img = highlightMolecules(mols, mcs1, 3) # SVGデータを取得 molsvg = img.data # 背景を透明に設定 molsvg = molsvg.replace("opacity:1.0", "opacity:0.0"); # ラベルのサイズを大きくする molsvg = molsvg.replace("12px", "20px") # 変更したSVGデータを保存する f = open("../data/T6/mcs_largestcluster.svg", "w") f.write(molsvg) f.close()
閾値を設定
部分構造の閾値を下げることも可能で、例えば、MCSが入力構造の80%の化合物にだけあれば良いようにすることもできます。
mcs2 = rdFMCS.FindMCS(mols, threshold=0.8) print('SMARTS string: %s'%mcs2.smartsString) # 部分構造を描画 m2 = Chem.MolFromSmarts(mcs2.smartsString) Draw.MolsToGridImage([m1,m2], legends=['mcs1', 'mcs2: thr 0.8']) # SMARTS string: [#6]=[#6]-[#6](=[#8])-[#7]-[#6]1:[#6]:[#6]2:[#6](-[#7]-[#6]3:[#6]:[#6]:[#6]:[#6]:[#6]:3-[#9]):[#7]:[#6]:[#7]:[#6]:2:[#6]:[#6]:1-[#8]-[#6]-[#6]-[#7]-[#6]
highlightMolecules(mols, mcs2, 12)
この例を見ればわかるように、化合物CHEMBL27685は除外され(閾値 0.8)、その結果、より大きな共通部分構造として、メタ位にフッ素置換したベンゼンとより短いアルキル鎖をもつ部分構造が見つかっています。
訳注(2020/04)
今回の実行結果ではオリジナルのノートと異なりCHEMBL27685は上図に含まれません。上図中ではCHEMBL3622640にはMCSが含まれていません。閾値0.8のMCS(mcs2)はエーテル側鎖の先のN原子までを含んでいますが、CHEMBL3622640ではN原子ではなくO原子となっています。MCSの閾値を下げることで、全ての分子には含まれないもののより大きな部分構造をMCSとして見出していることが確認できます。
訳注ここまで
環結合(ring bonds)のマッチング
このパラメーターが True
にセットすると、環の結合だけが他の環の結合とマッチングします。
mcs3 = rdFMCS.FindMCS(mols, ringMatchesRingOnly=True, threshold=0.8) # 部分構造を描画 m3 = Chem.MolFromSmarts(mcs3.smartsString) Draw.MolsToGridImage([m1,m2,m3], legends=['mcs1', 'mcs2: +thr 0.8','mcs3: +ringmatch'])
ここでは選択した閾値とパラメーターに依存して、少し異なるMCSが得られています。 RDKitのFMCSモジュールにはもっと多くのパラメーターのオプションが利用可能で、例えば考慮する原子、結合、価数のマッチングといったものがある、ということを注記しておきます。
さあ、より多様なセットをさっと見てみましょう:ChEMBLからダウンロードしたEGFR化合物群
データを高活性な化合物(pIC50>9)だけに限定し、このサブセットから最大共通骨格を発見しましょう。
# 全EGFRデータの読み込み mol_df=pd.read_csv('../data/T1/EGFR_compounds.csv', index_col=0) print(mol_df.shape) # pIC50 > 9 (IC50 > 1nM)の化合物のみを保持 mol_df=mol_df[mol_df.pIC50>9] print(mol_df.shape) # データフレームに分子を追加 PandasTools.AddMoleculeColumnToFrame(mol_df, 'smiles') mol_df.head(3) # (5428, 5) # (191, 5)
molecule_chembl_id | units | IC50 | smiles | pIC50 | ROMol | |
---|---|---|---|---|---|---|
0 | CHEMBL63786 | nM | 0.003 | Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 | 11.522879 | |
1 | CHEMBL53711 | nM | 0.006 | CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 | 11.221849 | |
2 | CHEMBL35820 | nM | 0.006 | CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC | 11.221849 |
# 選択した高活性化合物だけについて計算を実行します mols= [] for idx, row in mol_df.iterrows(): m = Chem.MolFromSmiles(row.smiles) m.SetProp("_Name",row.molecule_chembl_id) mols.append(m) # 時間の都合上、このセットから50化合物をランダムに取り出します rand_mols = [mols[i] for i in numpy.random.choice(range(len(mols)), size=50)] mcs1 = rdFMCS.FindMCS(rand_mols) print('SMARTS string1: %s'%mcs1.smartsString) mcs2 = rdFMCS.FindMCS(rand_mols, threshold=0.8) print('SMARTS string2: %s'%mcs2.smartsString) mcs3 = rdFMCS.FindMCS(rand_mols, ringMatchesRingOnly=True, threshold=0.8) print('SMARTS string3: %s'%mcs3.smartsString) # 部分構造を描画 m1 = Chem.MolFromSmarts(mcs1.smartsString) m2 = Chem.MolFromSmarts(mcs2.smartsString) m3 = Chem.MolFromSmarts(mcs3.smartsString) Draw.MolsToGridImage([m1,m2,m3], legends=['mcs1', 'mcs2: +thr 0.8','mcs3: +ringmatch']) # SMARTS string1: [#6](:,-[#6]:,-[#6]:,-[#6]):,-[#6]:,-[#6]:,-[#7]:,-[#6] # SMARTS string2: [#6]1:[#6]:[#6]:[#6]:[#6](:[#6]:1):,-[#7]:,-[#6]:[#7]:[#6]-,:[#7]-,:[#6](:[#6]:[#6]:[#6]):[#6] # SMARTS string3: [#6]:[#6]:[#6](:[#7]:[#6]:[#7]:[#6]-[#7]-[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1):[#6]
highlightMolecules(rand_mols, mcs2, 12)
from ipywidgets import interactive, interact, fixed, Text, widgets from IPython.display import SVG,display
def renderMCS(perc): tmcs = rdFMCS.FindMCS(rand_mols, threshold=perc/100.) if tmcs is None: print('No MCS found') return None else: m = Chem.MolFromSmarts(tmcs.smartsString) print(tmcs.smartsString) return(display(SVG(IPythonConsole._toSVG(m))))
# スライダーが反応するのに数秒かかるかもしれません w = interact(renderMCS, perc=widgets.IntSlider(min=0,max=100,step=20,value=80)) # interactive(children=(IntSlider(value=80, description='perc', step=20), Output()), _dom_classes=('widget-inter…
訳注(2020/04)
Jupyter Notebook場ではスライダーと画像が表示され、スライダーにより閾値を変更すると合わせてMCSが変化し画像が表示される、という構成になっています。Markdown形式には移行できなかったので興味を持たれた方はご自身の環境で試されてください。
訳注ここまで
クイズ
- 最大共通部分構造(MCS)の計算が役立つのはなぜでしょう?
- どうすればMCSを計算できるか簡潔に説明できますか?
- EGFR活性化合物の典型的なフラグメントはどのようにみえますか?
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。
MCSに関して:
* 化学の新しいカタチ:RDKitを用いた部分構造検索とMCSアルゴリズム
* AI創薬のためのケモインフォマティクス 入門:7章グラフ構造を利用した類似性の評価
化学構造とグラフ??という場合には以下の記事をお勧めします。
* 化学構造情報とグラフ理論
勝手にいろいろと引用させていただいておりますので問題があれば削除します。
TeachOpenCADD トピック5 〜化合物クラスタリング〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル5をもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下、訳。
トークトリアル 5
化合物クラスタリング
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Calvinna Caswara and Gizem Spriewald
このトークトリアルの目的
類似の化合物は同じターゲット分子に結合し、類似の効果を示すかもしれません。この類似性質原則(similar property principle、SPP)に基づき、化合物類似度をつかって、クラスタリングによる化合物集団(chemical group)の構築を行うことができます。そのようなクラスタリングによって、より大きなスクリーニング化合物セットから、さらなる実験を行う対象となる化合物の多様性のあるセットを選ぶことができます。
学習の目標
このトークトリアルでは、以下についてさらに学びます:
理論
実践
- Butinaクラスタリングと化合物選択の例
レファレンス
- Butina, D. Unsupervised Data Base Clustering Based on Daylight’s Fingerprint and Tanimoto Similarity: A Fast and Automated Way To Cluster Small and Large Data Set. J. Chem. Inf. Comput. Sci. 1999.
- Leach, Andrew R., Gillet, Valerie J. An Introduction to Chemoinformatics. 2003
- Jarvis-Patrickクラスタリング: http://www.improvedoutcomes.com/docs/WebSiteDocs/Clustering/Jarvis-Patrick_Clustering_Overview.htm
- TDT チュートリアル: https://github.com/sriniker/TDT-tutorial-2014/blob/master/TDT_challenge_tutorial.ipynb
- rdkitクラスタリングドキュメンテーション: http://rdkit.org/docs/Cookbook.html#clustering-molecules
理論
クラスタリングとJarvis-Patrickアルゴリズムの紹介
クラスタリング(Clustering) は「ある一連のものをグループ化するタスクで、(クラスターと呼ばれる)同じグループに属するものが(ある基準において)、他のグループ(クラスター)に属するものと比較して、お互いにより似通っているようにグループ化するタスク」と定義することができます。
医薬品研究における化合物クラスタリングは化合物間の化学的、構造的類似度にしばしば基づいており、共通する特徴をもつグループをみつけたり、さらに解析を行うための多様性がある代表的な化合物セットをデザインするために行われます。
一般的な手順:
- 近接するポイント間の類似度によってデータをクラスタリングすることに基づく手法
- ケモインフォマティクスにおいて、化合物は化合物フィンガープリントでエンコードされ、類似度はタニモト類似度によって表現されることが多いです(トークトリアル4参照)
- 念のため再確認:フィンガープリントはバイナリのベクトルで、各ビットはある特定の部分構造フラグメントが分子の中にあるか無いかを示します。
- 類似度(あるいは距離)行列: バイナリのフィンガープリントによって表現された分子のペアの類似度はタニモト係数を使って定量化されることがもっとも多く、これは共通する特徴(ビット)の数を評価する値です。
- タニモト係数の値の範囲は0(類似性なし)から1(類似度が高い)となっています。
利用可能なクラスタリングアルゴリズムはたくさんありますが、中でもJarvis-Patrickクラスタリング は医薬品の分野において最も広く使われているアルゴリズムの一つです。
Jarvis-Patrickクラスタリングアルゴリズムは2つのパラメーター、KとKminで定義されています:
- 各分子のK最近傍のセットを計算
- 2つの化合物が次の場合、同じクラスターに含める
- 最近傍のリストにお互い含まれている場合
- K最近傍のうち少なくともKminを共通に持つ場合
Jarvis-Patrickクラスタリングアルゴリズムは決定的アルゴリズムで、巨大な化合物セットを数時間のうちに処理することができます。しかしながら、この手法の悪い側面として、大きく不均質なクラスタを生み出す傾向があります(Butinaクラスタリングの参考文献を参照)。
他のクラスタリングアルゴリズムについて知りたければ scikit-learn clustering module を参照してください。
Butinaクラスタリングの詳細な説明
Butinaクラスタリング(J. Chem. Inf. Model.(1999), 39(4), 747)はより小さいが均質なクラスターを見つけるために開発されました。(少なくとも)クラスターの中心が、クラスター内の他の全ての分子と、与えられた閾値以上に類似していること、を前提条件とします。
このクラスタリングアプローチには鍵となるステップがあります(下のフローチャート参照してください):
1. データの準備と化合物エンコーディング
- 化合物類似性を見つけるために、入力データの化合物(例 SMILESとしてあたえられたもの)は化合物フィンガープリントとしてエンコードされます。例えば、RDK5フィンガープリントは、よく知られているDaylight Fingerprint(のオリジナル文献で使われていたもの)に似たサブグラフに基づく(subgraph-based)フィンガープリントです。
2. タニモト類似度(あるいは距離)行列
- 2つのフィンガープリント間の類似度はタニモト係数を使って計算します
- 全ての可能な分子/フィンガープリントのペア間のタニモト類似度を要素とする行列です( n x n 類似度行列、n=化合物数、上三角行列のみが使われます)
- 同様に、距離行列が計算できます(1 - 類似度)
3. 化合物クラスタリング:中心(centroids)と排除球(exclusion spheres)
注:化合物は次の場合に同じクラスターに含められます。(距離行列を使用する場合)クラスターのセントロイドからの最長距離が特定のカットオフ値以下の時、あるいは、(類似度行列を使用する場合)最小の類似度が特定のカットオフ値以上のときです。
排除球(exclusion shpheres)に基づくクラスタリング
- 並べ替えたリスト内の最初の分子(中心、centroid)から始めます。
- クラスタリング手順の最後までフラグの立てられなかった化合物はシングルトン(singleton)となります。
- シングルトンとして割り当てられた化合物の中には、指定したタニモト類似度インデックスで近傍となる化合物をもつものありえますが、これらの近傍化合物は、クラスター中心に基づく枠組みによってすでに排除球に含められていることに注意してください。
from IPython.display import IFrame IFrame('./images/butina_full.pdf', width=600, height=300)
Figure 1: Butinaクラスタリングアルゴリズムの理論的な例(Calvinna Caswaraによる)
多様な化合物の選択
代表的な化合物セットを見つけるという考え方が、医薬品業界ではしばしば使われます。
- 例えば、バーチャルスクリーニングキャンペーンを行なったものの、リソースが非常に限られており、確認のためのアッセイ実験で評価可能な化合物数が少ないとしましょう。
- このスクリーニングから得られる情報を最大限のものとするために、多様性のある化合物セットを選びたいと思います。そのためには、活性を有する可能性がある化合物のリストの中で、各ケミカルシリーズを代表する一つを選択します。
別のシナリオとしては、構造活性相関、即ち化合物の小さな構造変化がin vitro活性にどう影響を与えるか?、について情報を得るために一連の化合物を選択する、というシナリオが考えられます。
実践
Butinaクラスタリングアルゴリズムの使用例
S. Riniker と G. LandrumによるTDT tutorial notebookの例にしたがって適用してみます。
1. データの読み込みとフィンガープリントの計算
ここではデータの準備とフィンガープリントの計算を行います。
# パッケージのimport import pandas as pd import numpy import matplotlib.pyplot as plt import time import random from random import choices from rdkit import Chem from rdkit.Chem import AllChem from rdkit import DataStructs from rdkit.DataStructs import cDataStructs from rdkit.ML.Cluster import Butina from rdkit.Chem import Draw from rdkit.Chem import rdFingerprintGenerator from rdkit.Chem.Draw import IPythonConsole
# データを読み込んで眺める # トークトリアル2で取り出したフィルタリング済みのデータ compound_df= pd.read_csv('../data/T2/EGFR_compounds_lipinski.csv',sep=";", index_col=0) print('data frame shape:',compound_df.shape) compound_df.head() # data frame shape: (4523, 10)
molecule_chembl_id | units | IC50 | smiles | pIC50 | MW | HBA | HBD | LogP | rule_of_five_conform | |
---|---|---|---|---|---|---|---|---|---|---|
0 | CHEMBL63786 | nM | 0.003 | Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 | 11.522879 | 349.021459 | 3 | 1 | 5.2891 | yes |
1 | CHEMBL53711 | nM | 0.006 | CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 | 11.221849 | 343.043258 | 5 | 1 | 3.5969 | yes |
2 | CHEMBL35820 | nM | 0.006 | CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC | 11.221849 | 387.058239 | 5 | 1 | 4.9333 | yes |
3 | CHEMBL53753 | nM | 0.008 | CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 | 11.096910 | 329.027607 | 5 | 2 | 3.5726 | yes |
4 | CHEMBL66031 | nM | 0.008 | Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 | 11.096910 | 339.011957 | 4 | 2 | 4.0122 | yes |
# SMILESから化合物を作成し、配列に格納 mols = [] for i in compound_df.index: chemblId = compound_df['molecule_chembl_id'][i] smiles = compound_df['smiles'][i] mols.append((Chem.MolFromSmiles(smiles), chemblId)) mols[0:5] # [(<rdkit.Chem.rdchem.Mol at 0x1037674e0>, 'CHEMBL63786'), # (<rdkit.Chem.rdchem.Mol at 0x103767bc0>, 'CHEMBL53711'), # (<rdkit.Chem.rdchem.Mol at 0x103767df0>, 'CHEMBL35820'), # (<rdkit.Chem.rdchem.Mol at 0x103767da0>, 'CHEMBL53753'), # (<rdkit.Chem.rdchem.Mol at 0x103767b70>, 'CHEMBL66031')]
# 全化合物についてフィンガープリントを作成 rdkit_gen = rdFingerprintGenerator.GetRDKitFPGenerator(maxPath=5) fingerprints = [rdkit_gen.GetFingerprint(m) for m,idx in mols] # 化合物/フィンガープリントは幾つあるでしょうか? print('Number of compounds converted:',len(fingerprints)) print('Fingerprint length per compound:',len(fingerprints[0])) # Number of compounds converted: 4523 # Fingerprint length per compound: 2048
2. タニモト類似度と距離行列
フィンガープリントが生成できたので、次の段階に進みます。クラスター中心となる可能性があるものの同定です。このために、タニモト類似度と距離行列を計算する関数を定義します。
# フィンガープリントのリストのための距離行列を計算 def Tanimoto_distance_matrix(fp_list): dissimilarity_matrix = [] for i in range(1,len(fp_list)): similarities = DataStructs.BulkTanimotoSimilarity(fp_list[i],fp_list[:i]) # 距離行列が必要なので、類似度行列の全要素について1-xを計算 dissimilarity_matrix.extend([1-x for x in similarities]) return dissimilarity_matrix
RDKitクックブック:分子のクラスタリングも参照してください。
# 例:2つのフィンガープリントの類似度を一つ計算 sim = DataStructs.TanimotoSimilarity(fingerprints[0],fingerprints[1]) print ('Tanimoto similarity: %4.2f, distance: %4.2f' %(sim,1-sim)) # Tanimoto similarity: 0.68, distance: 0.32
# 例:距離行列を計算(距離 = 1-類似度) Tanimoto_distance_matrix(fingerprints)[0:5] # [0.31818181818181823, # 0.31474103585657376, # 0.430976430976431, # 0.26991150442477874, # 0.07111111111111112]
# 追記: 行列ではなくリストのように見えます。 # しかし、リスト形式の上三角行列です。 n = len(fingerprints) # n*(n-1)/2 によって三角行列の要素数を計算します。 elem_triangular_matr = (n*(n-1))/2 print(int(elem_triangular_matr), len(Tanimoto_distance_matrix(fingerprints))) # 10226503 10226503
3. 化合物クラスタリング:中心(centroids)と排除球(exclusion spheres)
ここでは、化合物のクラスター化を行い結果を眺めます。
クラスタリング関数の定義
# 入力:フィンガープリントとクラスタリングの閾値 def ClusterFps(fps,cutoff=0.2): # タニモト距離行列の計算 distance_matr = Tanimoto_distance_matrix(fps) # Butinaアルゴリズムの実装を使ってデータをクラスタリング: clusters = Butina.ClusterData(distance_matr,len(fps),cutoff,isDistData=True) return clusters
フィンガープリントの類似度に基づき化合物をクラスター化
# データセットのためのクラスタリング手順を実行 clusters = ClusterFps(fingerprints,cutoff=0.3) # クラスターの数とサイズについて短いレポートを表示 num_clust_g1 = len([c for c in clusters if len(c) == 1]) num_clust_g5 = len([c for c in clusters if len(c) > 5]) num_clust_g25 = len([c for c in clusters if len(c) > 25]) num_clust_g100 = len([c for c in clusters if len(c) > 100]) print("total # clusters: ", len(clusters)) print("# clusters with only 1 compound: ", num_clust_g1) print("# clusters with >5 compounds: ", num_clust_g5) print("# clusters with >25 compounds: ", num_clust_g25) print("# clusters with >100 compounds: ", num_clust_g100) # total # clusters: 681 # # clusters with only 1 compound: 305 # # clusters with >5 compounds: 157 # # clusters with >25 compounds: 33 # # clusters with >100 compounds: 4
訳注(04/2020)
ここで実行しているRDKitのButinaモジュールの戻り値はクラスターの情報を含むタプルのタプルです。
( (cluster1_elem1, cluster1_elem2, …), (cluster2_elem1, cluster2_elem2, …), …)
という形式で、タプル内タプルの最初のclusterX_elem1
は各クラスターXのcentroidとなっています。 上記ではlen()
でクラスターのサイズを取得しています。
訳注ここまで
# クラスターの大きさをプロット fig = plt.figure(1, figsize=(10, 4)) plt1 = plt.subplot(111) plt.axis([0, len(clusters), 0, len(clusters[0])+1]) plt.xlabel('Cluster index', fontsize=20) plt.ylabel('Number of molecules', fontsize=20) plt.tick_params(labelsize=16) plt1.bar(range(1, len(clusters)), [len(c) for c in clusters[:len(clusters)-1]], lw=0) plt.show()
合理的なカットオフ値を選ぶにはどうすればよいか?
クラスタリングの結果はユーザーの選んだ閾値に依存するので、カットオフ値の選択についてより詳細に見てみたいと思います。
for i in numpy.arange(0., 1.0, 0.1): clusters = ClusterFps(fingerprints,cutoff=i) fig = plt.figure(1, figsize=(10, 4)) plt1 = plt.subplot(111) plt.axis([0, len(clusters), 0, len(clusters[0])+1]) plt.xlabel('Cluster index', fontsize=20) plt.ylabel('Number of molecules', fontsize=20) plt.tick_params(labelsize=16) plt.title('Threshold: '+str('%3.1f' %i), fontsize=20) plt1.bar(range(1, len(clusters)), [len(c) for c in clusters[:len(clusters)-1]], lw=0) plt.show()
結果を見ると、閾値(距離のカットオフ)が高いほど、より多くの化合物が類似しているとみなされ、したがってより少ないクラスター数にクラスタリングされます。閾値が低くなると、より多数の小さなクラスターとシングルトン(singleton)が現れます。
- 距離のカットオフ値がより小さくなるほど、一つのクラスターに所属する化合物は互いにより類似していることが求められます。
上のプロットをもとにて、私たちは距離の閾値0.2を選択しました。シングルトンの数が多くなく、クラスターのサイズが大き過ぎていませんが、分布は滑らかになっています。
dist_co = 0.2 clusters = ClusterFps(fingerprints,cutoff=dist_co) # クラスターのサイズをプロットー保存 fig = plt.figure(1, figsize=(8, 2.5)) plt1 = plt.subplot(111) plt.axis([0, len(clusters), 0, len(clusters[0])+1]) plt.xlabel('Cluster index', fontsize=20) plt.ylabel('# molecules', fontsize=20) plt.tick_params(labelsize=16) plt1.bar(range(1, len(clusters)), [len(c) for c in clusters[:len(clusters)-1]], lw=0) plt.title('Threshold: '+str('%3.1f' %dist_co), fontsize=20) plt.savefig("../data/T5/cluster_dist_cutoff_%4.2f.png" %dist_co, dpi=300, bbox_inches="tight", transparent=True) print('Number of clusters %d from %d molecules at distance cut-off %4.2f' %(len(clusters), len(mols), dist_co)) print('Number of molecules in largest cluster:', len(clusters[0])) print('Similarity between two random points in same cluster %4.2f'%DataStructs.TanimotoSimilarity(fingerprints[clusters[0][0]],fingerprints[clusters[0][1]])) print('Similarity between two random points in different cluster %4.2f'%DataStructs.TanimotoSimilarity(fingerprints[clusters[0][0]],fingerprints[clusters[1][0]])) # Number of clusters 1132 from 4523 molecules at distance cut-off 0.20 # Number of molecules in largest cluster: 143 # Similarity between two random points in same cluster 0.92 # Similarity between two random points in different cluster 0.30
クラスターの可視化
最大のクラスターからの化合物例10個
最初の最もサイズの大きいクラスターから、初めの10化合物をより詳細に見てみましょう。
print ('Ten molecules from largest cluster:') # 化合物描画 Draw.MolsToGridImage([mols[i][0] for i in clusters[0][:10]], legends=[mols[i][1] for i in clusters[0][:10]], molsPerRow=5) # Ten molecules from largest cluster:
# トークトリアル9でMCS解析を行うため最大のサイズのクラスターの化合物を保存 # 訳注:上記はトークトリアル6の記載ミスか? w = Chem.SDWriter('../data/T5/molSet_largestCluster.sdf') # データの準備 tmp_mols=[] for i in clusters[0]: tmp = mols[i][0] tmp.SetProp("_Name",mols[i][1]) tmp_mols.append(tmp) # データの書き出し for m in tmp_mols: w.write(m)
2番目に大きいクラスターからの化合物例10個
print ('Ten molecules from second largest cluster:') # 化合物の描画 Draw.MolsToGridImage([mols[i][0] for i in clusters[1][:10]], legends=[mols[i][1] for i in clusters[1][:10]], molsPerRow=5) # Ten molecules from second largest cluster:
それぞれのクラスターの最初の10化合物は実際にお互いに類似しており、多くが共通の骨格を共有しています(視覚的に判断する限り)。
化合物セットの最大共通部分構造(maximum common substructure、MCS)を計算する方法について、もっと知りたければトークトリアル6を参照してください。
最初の10クラスターからの例
比較のため、最初の10クラスターのクラスター中心を見てみます。
print ('Ten molecules from first 10 clusters:') # 化合物を描画 Draw.MolsToGridImage([mols[clusters[i][0]][0] for i in range(10)], legends=[mols[clusters[i][0]][1] for i in range(10)], molsPerRow=5) # Ten molecules from first 10 clusters:
最初の3つクラスターから取り出したクラスター中心をSVGファイルとして保存。
# イメージを生成 img = Draw.MolsToGridImage([mols[clusters[i][0]][0] for i in range(0,3)], legends=["Cluster "+str(i) for i in range(1,4)], subImgSize=(200,200), useSVG=True) # SVGデータを取得 molsvg = img.data # 不透明な背景を透明に置き換え、フォントサイズを設定 molsvg = molsvg.replace("opacity:1.0", "opacity:0.0"); molsvg = molsvg.replace("12px", "20px"); # 変換したSVGデータをファイルに保存 f = open("../data/T5/cluster_representatives.svg", "w") f.write(molsvg) f.close()
まだいくらか類似性が残っているのが見てとれますが、明らかに、異なるクラスターの中心同士は、一つのクラスター内の化合物同士と比較して、類似性が低くなっています。
クラスター内タニモト類似度
クラスター内のタニモト類似度をみてみることもできます。
# 各クラスターのフィンガープリントの全てのペアについてタニモト類似度を計算する関数 def IntraTanimoto(fps_clusters): intra_similarity =[] # クラスターごとの内部類似度を計算 for k in range(0,len(fps_clusters)): # 類似度行列(1ー距離)に変換されたタニモト距離行列の関数 intra_similarity.append([1-x for x in Tanimoto_distance_matrix(fps_clusters[k])]) return intra_similarity
# 最初の10個のクラスターについてフィンガープリントを再計算 mol_fps_per_cluster=[] for c in clusters[:10]: mol_fps_per_cluster.append([rdkit_gen.GetFingerprint(mols[i][0]) for i in c]) # クラスター内類似度を計算 intra_sim = IntraTanimoto(mol_fps_per_cluster)
# クラスター内類似度とともにバイオリンプロット pos = list(range(10)) labels = pos plt.figure(1, figsize=(10, 5)) ax = plt.subplot(111) r = plt.violinplot(intra_sim, pos, showmeans=True, showmedians=True, showextrema=False) ax.set_xticks(pos) ax.set_xticklabels(labels) ax.set_yticks(numpy.arange(0.6, 1., 0.1)) ax.set_title('Intra-cluster Tanimoto similarity', fontsize=13) r['cmeans'].set_color('red') # 平均=赤色, 中央値=青色
化合物選択
以下では、多様性のあるサブセットとして最大 1000化合物の最終リストを取り出します。
このため、最大サイズのクラスターから始めて最大1000化合物を取り出すまで、各クラスターのクラスター中心(即ち、各クラスターの最初の分子)をとりだして、各クラスターのクラスター中心と最も似ている10化合物(あるいは、クラスター化合物が10以下のときは50%の化合物)を選びます。これにより、各クラスターの代表が得られます。
この化合物選択の目的は、確認用アッセイ実験に提案するためのより小さな化合物セットの多様性を確保することです。
選択の手順はS. Riniker と G. LandrumによるTDT tutorial notebookからとりました。ノートブックで述べられているように、このアプローチの背景にあるアイデアは、確認アッセイ実験の結果からSARを取得しつつ(より大きなクラスターの非常に類似した化合物群を保ちつつ)、多様性を確保する(各クラスターの代表を組み入れる)ことです。
クラスター中心を取得します。
# 各クラスターのクラスター中心を取得(各クラスターの最初の化合物) clus_center = [mols[c[0]] for c in clusters] # クラスター中心/クラスターはいくつあるか? print('Number of cluster centers: ', len(clus_center)) # Number of cluster centers: 1132
サイズによってクラスターを並べ替え、類似度によって各クラスターの化合物を並べ換える。
# クラスター中心に対する類似度に基づいてクラスター内の化合物を並べ替え、 # サイズに基づきクラスターを並べ替えます。 clusters_sort = [] for c in clusters: if len(c) < 2: continue # シングルトン else: # 各クラスター要素についてフィンガープリントを計算 fps_clust = [rdkit_gen.GetFingerprint(mols[i][0]) for i in c] # クラスター中心に対する全クラスターメンバーの類似度 simils = DataStructs.BulkTanimotoSimilarity(fps_clust[0],fps_clust[1:]) # (中心は除いて!)類似度に化合物のインデックスを付与する simils = [(s,index) for s,index in zip(simils, c[1:])] # 類似度によって降順に並べ替え simils.sort(reverse=True) # クラスターのサイズと化合物インデックスをclusters_sortに保存 clusters_sort.append((len(simils), [i for s,i in simils])) # クラスターサイズによって降順にソート clusters_sort.sort(reverse=True)
最大1000化合物を取り出す。
# 選んだ分子を数える、はじめにクラスター中心を取り出す sel_molecules = clus_center.copy() # 最大のクラスターからはじめて、各クラスターの10化合物(あるいは 最大の50%)を取り出す index = 0 diff = 1000 - len(sel_molecules) while diff > 0 and index < len(clusters_sort): # 並べ替えたクラスターのインデックスを取得 tmp_cluster = clusters_sort[index][1] # 最初のクラスターが10以上の大きさの時、ちょうど10化合物を取り出す if clusters_sort[index][0] > 10: num_compounds = 10 # 10より小さいなら、化合物の半分をとる else: num_compounds = int(0.5*len(c))+1 if num_compounds > diff: num_compounds = diff # picked_fps と名付けたリストのリストに取り出した分子と構造を書き込む # 訳注: 上記リストの名称 picked_fpsはsel_moleculesの記載ミスか? sel_molecules += [mols[i] for i in tmp_cluster[:num_compounds]] index += 1 diff = 1000 - len(sel_molecules) print('# Selected molecules: '+str(len(sel_molecules))) # # Selected molecules: 1132
この多様性のある化合物のセットは実験による評価に用いることができます。
訳注(2020/04)
上記で選択された化合物は1132個となっており、意図した最大1000個とはなっていません。オリジナルのノートブックと異なり、現在(2020年04月)の取得データをもとにすると、クラスター数が1132個となり、クラスター数の時点で1000個を超えています(オリジナルは988個)。ですので、上記コードのクラスター中心を取り出したsel_molecules
が1000以上(diff <0
) となり、while
条件文以降が実行されないという結果となっています。以下、コードが正しく機能するか確認のため、取り出す化合物の最大数を1200個にして実行して見ます。
sel_molecules2 = clus_center.copy() index = 0 diff = 1200 - len(sel_molecules2) while diff > 0 and index < len(clusters_sort): tmp_cluster = clusters_sort[index][1] if clusters_sort[index][0] > 10: num_compounds = 10 else: num_compounds = int(0.5*len(c))+1 if num_compounds > diff: num_compounds = diff sel_molecules2 += [mols[i] for i in tmp_cluster[:num_compounds]] index += 1 diff = 1200 - len(sel_molecules2) print('# Selected molecules: '+str(len(sel_molecules2))) # # Selected molecules: 1200
指定した1200個の化合物を取り出すことができました。他の検証方法としてはクラスタリングの閾値を0.2よりも大きくしてクラスター数を1000個以下とすることがあげられると思います。
訳注ここまで
(追加情報:実行時間)
トークトリアルの終わりに、データセットのサイズが変わった時に、Butinaクラスタリングの実行時間がどのように変わるかみてみましょう。
# 古いデータセットを再利用
sampled_mols = mols.copy()
より大きなデータセットを試してみることもできますが、10000以上のデータポイントの時点で非常に大きなメモリと時間がかかりはじめます(これがここで止めた理由です。)
# 時間を図るためのヘルパー関数 def MeasureRuntime(sampled_mols): start_time = time.time() sampled_fingerprints = [rdkit_gen.GetFingerprint(m) for m,idx in sampled_mols] # データセットでクラスタリングを実行 sampled_clusters = ClusterFps(sampled_fingerprints,cutoff=0.3) return(time.time() - start_time)
dsize=[100, 500, 1000, 2000, 4000, 6000, 8000, 10000] runtimes=[] # 置き換えありで、ランダムにサンプルをとる for s in dsize: tmp_set = [sampled_mols[i] for i in sorted(numpy.random.choice(range(len(sampled_mols)), size=s))] tmp_t= MeasureRuntime(tmp_set) print('Dataset size %d, time %4.2f seconds' %(s, tmp_t)) runtimes.append(tmp_t) # Dataset size 100, time 0.10 seconds # Dataset size 500, time 0.49 seconds # Dataset size 1000, time 0.99 seconds # Dataset size 2000, time 2.31 seconds # Dataset size 4000, time 5.96 seconds # Dataset size 6000, time 11.26 seconds # Dataset size 8000, time 18.28 seconds # Dataset size 10000, time 25.67 seconds
plt.plot(dsize, runtimes, 'g^') plt.title('Runtime measurement of Butina Clustering with different dataset sizes') plt.xlabel('# Molecules in data set') plt.ylabel('Runtime in seconds') plt.show()
クイズ
- 化合物のクラスタリングはなぜ重要なのでしょうか?
- 化合物セットのクラスタリングにどのアルゴリズムを使うことができますか?また、アルゴリズムの背景になる一般的な考え方はどのようなものでしょうか?
- 他のクラスタリングアルゴリズムを知っていますか?
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。
kNNによるクラスタリング:
* 化学の新しいカタチ:RDKitとscikit-learnで機械学習:変異原性をk-最近傍法で予測
Butinaクラスタリングは以下でも取り上げられています:
* AI創薬のためのケモインフォマティクス 入門:6章化合物の類似性を評価してみる
多様なライブラリの構築とその方法については以下で説明してくださっています:
* 化学の新しいカタチ:ケモインフォマティクスで多様な化合物ライブラリーを構築する
また、今回のトークトリアルでは取り上げられていませんでしたが、クラスタリングと関連する可視化として次元圧縮によるプロットがあげられると思います。ケミカルスペースの可視化については以下で解説してくださっています。
勝手にいろいろと引用させていただいておりますので問題があれば削除します。
TeachOpenCADD トピック4 〜リガンドベーススクリーニング〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル4をもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!
- トークトリアル 4
- リガンドベーススクリーニング:化合物類似性
- このトークトリアルの目的
- 学習の目標
- レファレンス
- 理論
- 実践
- ディスカッション
- クイズ
- 追記
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下、訳。
トークトリアル 4
リガンドベーススクリーニング:化合物類似性
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Andrea Morger and Franziska Fritz
このトークトリアルの目的
このトークトリアルでは、化合物をエンコード(記述子、フィンガープリント)し、比較(類似性評価)する様々なアプローチを取り扱います。さらに、バーチャルスクリーニングを実施します。バーチャルスクリーニングは、ChEMBLデータベースから取得し、リピンスキーのルールオブファイブでフィルタリングをかけた、EGFRに対して評価済みの化合物データセット(トークトリアル 2 参照) に対して、EGFR阻害剤ゲフィチニブ(Gefitinib)との類似性検索を実施するという形で実施します。
学習の目標
理論
- 化合物類似性(Molecular similarity)
- 化合物記述子(Molecular descriptors)
- 化合物フィンガープリント(Molecular fingerprints)
- 部分構造ベースのフィンガープリント(Substructure-based fingerprints)
- MACCSフィンガープリント(MACCS fingerprints)
- Morganフィンガープリント、サーキュラーフィンガープリント(Morgan fingerprints, circular fingerprints)
- 化合物類似性評価(Molecular similarity measures)
- タニモト係数(Tanimoto coefficient)
- Dice係数(Dice coefficient)
- バーチャルスクリーニング(Virtual screening)
- 類似性検索(similarity search)によるバーチャルスクリーニング
実践
- 分子の読み込みと描画
- 化合物記述子の計算
- 1D 化合物記述子:分子量
- 2D 化合物記述子:MACCS ファインガープリント
- 2D 化合物記述子:Morgan フィンガープリント
- 化合物類似性の計算
- MACCS フィンガープリント:タニモト類似性とDice類似性
- Morgan フィンガープリント:タニモト類似性とDice類似性
- 類似性検索によるバーチャルスクリーニング
- データセットの全化合物に対する化合物クエリの比較
- 類似度の分布
- 最も類似した分子の描画
- エンリッチメントプロットの作成
レファレンス
- レビュー"Molecular similarity in medicinal chemistry" (J. Med. Chem. (2014), 57, 3186-3204)
- RDKitのMorganフィンガープリント (RDKit tutorial on Morgan fingerprints)
- ECFP - extended-connectivity fingerprints (J. Chem. Inf. Model. (2010), 50,742-754)
- ケミカルスペース (ACS Chem. Neurosci. (2012), 19, 649-57)
- RDKitの化合物記述子リスト (RDKit documentation: Descriptors)
- RDKitのフィンガープリントのリスト (RDKit documentation: Fingerprints)
- エンリッチメントプロット(Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 313-31)
理論
化合物類似性
化合物類似性は化学情報学(ケモインフォマティクス、chemical informatics)の中でよく知られており、頻繁に用いられる考え方です。化合物とその特性の比較はいろいろな用途で使用でき、望みの特性と生理活性をもつ新しい化合物を見つけるのに役立つかもしれません。
構造的に類似の化合物は類似の特性、そして類似の生理活性を示すという考え方は、類似性質原則(similar property principle、SPP)や構造活性相関(structure activity relationship、SAR)に表れています。この文脈において、バーチャルスクリーニングは、結合親和性のわかっている化合物セットがあれば、そのような化合物をさらに探すことができる、というアイデアに基づいています。
化合物記述子
類似度は適用範囲に応じて様々な方法で評価することができます(参考 J. Med. Chem. (2014), 57, 3186-3204):
- 1D 化合物記述子: 溶解度、logP、分子量、融点 etc.
- 2D 化合物記述子: 分子グラフ(Molecular graph)、経路(path)、フラグメント、原子環境(atom environment)
- 分子の個々の部位の詳細な表現
- 一つの分子に対して多数のフィンガープリントと呼ばれる特徴/ビット
- 類似性検索と機械学習で非常によく使われる
- 3D 化合物記述子: 形状(Shape), 立体化学
- 化学者は通常2次元表現で訓練されている
- 化合物の自由度(flexibility、化合物の「正しい」配座はどれか?)のため、2次元表現と比べて頑健性が低い
- 化学者は通常2次元表現で訓練されている
- 生物学的類似性
- 生物学的フィンガープリント(例、個々のビットが異なるターゲット分子に対して評価された生理活性を表す)
- 化合物構造からは独立
- 実験データ(あるいは予測値)が必要
すでに トークトリアル 2 で、分子量やlogPといった1D 物理化学パラメーターを計算する方法を学びました。RDKitに実装されているそのような記述子は RDKit documentation: Descriptors で見つけることができます。
以降では、2D(あるいは3D)化合物記述子の定義に焦点を当てます。多くの場合、分子ごとに固有の(ユニークな)ものとなるので、これらの記述子はフィンガープリント(指紋)ともよばれます。
化合物フィンガープリント(Molecular fingerprints)
部分構造に基づくフィンガープリント(Substructure-based fingerprints)
化合物フィンガープリントは化学的な特徴と分子の特徴をビット文字列(bitstring)やビットベクトル(bitvector)、配列(array)の形でエンコードします。各ビットは、事前に定義された分子の特徴あるいは環境に相当し、「1」はその特徴が存在していることを、「0」は存在していないことを示します。実装の仕方によっては、数え上げベース(count-based)となっていて、ある特定の特徴がいくつ存在しているかを数えるようになっていることに注意してください。
フィンガープリントのデザインには複数の方法があります。ここではよく使われる2Dフィンガープリントのとして、MACCSキーとMorganフィンガープリントの2種類を導入します。 RDKit documentation: Fingerprintsに記載されているように、RDKitではこの2つ以外にも多数のフィンガープリントを提供しています。
MACCS フィンガープリント(MACCS fingerprints)
Molecular ACCess System (MACCS) フィンガープリント、あるいはMACCS構造キーとも名付けられている手法は、あらかじめ定義された166個の構造フラグメントから構成されています。各位置は、ある特定の構造フラグメントあるいはキーが存在しているかいないかを問い合わせた(クエリ)結果を格納しています。それぞれのキーは創薬化学者によって経験的に定義されたもので、利用、解釈が容易です。(RDKit documentation: MACCS keys).
Figure 2: MACCSフィンガープリントの図例(Andrea Morgerによる図)
Morganフィンガープリントとサーキュラーフィンガープリント(Morgan fingerprints and circular fingerprints)
この一連のフィンガープリントはMorganアルゴリズムに基づいています。ビットは分子の各原子の円形状の環境(circular environment)に相当しています。半径(radius)によって、環境として考慮にいれる近接の結合と原子の数を設定します。ビット文字列の長さを定義することもでき、より長いビット文字列を希望する長さに縮められます。従って、Morganフィンガープリントはある特定の数のビットには制限されません。Morganフィンガープリントに関してもっと知りたい場合はRDKit documentation: Morgan fingerprints を参照してください。Extended connectivity fingerprints (ECFP)もよく使われるフィンガープリントで、Morganアルゴリズムのバリエーションから導かれています。さらなる情報は(J. Chem. Inf. Model. (2010), 50,742-754)を参照してください。
Figure 3: Morganサーキュラーフィンガープリントの図例(Andrea Morgerによる図)
化合物類似性評価
記述子/フィンガープリントの計算ができれば、それらを比較することで、二つの分子の間の類似度を評価することができます。化合物類似度は様々な類似度係数で定量化することができますが、よく使われる2つの指標はタニモト係数とDice係数です(Tanimoto and Dice index) (J. Med. Chem. (2014), 57, 3186-3204)。
タニモト係数(Tanimoto coefficient)
a: 化合物Aに存在する特徴の数
b: 化合物Bに存在する特徴の数
c: 化合物AとBで共有されている特徴の数
Dice係数(Dice coefficient)
a: 化合物Aに存在する特徴の数
b: 化合物Bに存在する特徴の数
c: 化合物AとBで共有されている特徴の数
類似度評価は通常、それぞれのフィンガープリントの正の(1の)ビットの数と、両者が共通してもつ正のビットの数を考慮します。Dice類似度は通常タニモト類似度よりも大きな値を返し、それはそれぞれの分母の違いに由来します。:
バーチャルスクリーニング(Virtual screening)
医薬品探索の初期段階における課題は、低分子(化合物)のセットを、有りうる巨大なケミカルスペースから、研究対象のターゲット分子に結合するポテンシャルのあるものに範囲を狭めることです。このケミカルスペースは非常に大きく、低分子化合物群は部分構造(chemical moiety)の1020 の組み合わせにまでいたります (ACS Chem. Neurosci. (2012), 19, 649-57) 。
目的のターゲット分子に対するこれら低分子の活性を評価するハイスループットスクリーニング(HTS)実験は費用と時間が非常にかかるので、計算機に支援された(computer-aided)手法により、試験にかける低分子のリストをより絞り込む(focused list)ことが期待されています。このプロセスはバーチャル(ハイスループット)スクリーニングと呼ばれていて、研究対象のターゲット分子に結合する可能性の最も高い低分子を見つけるために、巨大な低分子ライブラリーをルールとパターンのどちらか、あるいは両方によってフィルタリングします。
類似度検索を用いたバーチャルスクリーニング
バーチャルスクリーニングの簡単な方法として、既知の活性化合物(群)と新規な化合物セットを比較して、最も類似しているものを探すことが行われます。類似性質原則(similar property principle、SPP)に基づき、(例えば既知の阻害剤に)最も類似した化合物は類似の効果を有すると推測されます。類似性検索に必要となるものは次の通りです(上でも詳細に議論しています)。
- 化学/分子の特徴をエンコードした表現
- 特徴のポテンシャルの重み付け(オプション)
- 類似度評価
類似性検索はある特定のデータベースの全ての化合物と一つの化合物との間の類似度を計算することで実行することができます。データベースの化合物の類似度係数によるランク付けにより、最も類似度の高い分子が得られます。
エンリッチメントプロット(Enrichment plots)
エンリッチメントプロットはバーチャルスクリーニングの結果の妥当性を評価するために使われ、ランク付けされたリストの上位x%の中に含まれる活性化合物の比率を表します。すなわち、
[:350]
Figure 4: バーチャルスクリーニングの結果のエンリッチメントプロットの例
実践
実践編の最初のパートでは、RDKitを使って化合物のエンコード(化合物フィンガープリント)をしたのち、上の理論編で議論したように、類似度(化合物類似性評価)を計算するため、それらの比較を実行します。
2番目のパートではこれらのエンコーディングと比較手法を使って類似度検索(バーチャルスクリーニング)を実施します。既知のEGFR阻害剤ゲフィチニブ(Gefitinib)をクエリとして使用し、EGFRに対して試験済みの化合物データセットの中から類似した化合物を検索します。このデータセットは トークトリアル1でChEMBLから抽出し、トークトリアル2でリピンスキーのルールオブファイブによりフィルタリングをおこなったものです。
化合物の読み込みと描画
まず、8個の化合物例を定義し描画します。後ほど、これらの分子をエンコードし比較します。SMILES形式の分子をRDKitのmolオブジェクトに変換し、RDKitのDraw
関数で可視化します。
# 関連するPythonパッケージのimport # 基本的な分子を取り扱う機能はモジュール rdkti.Chem にあります from rdkit import Chem # 描画関係 from rdkit.Chem.Draw import IPythonConsole from rdkit.Chem import Draw from rdkit.Chem import Descriptors from rdkit.Chem import AllChem from rdkit.Chem import MACCSkeys from rdkit.Chem import rdFingerprintGenerator from rdkit import DataStructs import math import numpy as np import pandas as pd from rdkit.Chem import PandasTools import matplotlib.pyplot as plt
# SMILES形式の化合物 smiles1 = 'CC1C2C(C3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4O)O)O)O)C(=O)N)N(C)C)O' # Doxycycline smiles2 = 'CC1(C(N2C(S1)C(C2=O)NC(=O)C(C3=CC=C(C=C3)O)N)C(=O)O)C' # Amoxicilline smiles3 = 'C1=COC(=C1)CNC2=CC(=C(C=C2C(=O)O)S(=O)(=O)N)Cl' # Furosemide smiles4 = 'CCCCCCCCCCCC(=O)OCCOC(=O)CCCCCCCCCCC' # Glycol dilaurate smiles5 = 'C1NC2=CC(=C(C=C2S(=O)(=O)N1)S(=O)(=O)N)Cl' # Hydrochlorothiazide smiles6 = 'CC1=C(C(CCC1)(C)C)C=CC(=CC=CC(=CC(=O)O)C)C' # Isotretinoine smiles7 = 'CC1(C2CC3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4O)O)O)O)C(=O)N)N(C)C)O' # Tetracycline smiles8 = 'CC1C(CC(=O)C2=C1C=CC=C2O)C(=O)O' # Hemi-cycline D # 化合物SMILESのリストを作成 smiles = [smiles1, smiles2, smiles3, smiles4, smiles5, smiles6, smiles7, smiles8] # ROMolオブジェクトのリストを作成 mols = [Chem.MolFromSmiles(i) for i in smiles] # 化合物名称のリストを作成 mol_names = ['Doxycycline', 'Amoxicilline', 'Furosemide', 'Glycol dilaurate', 'Hydrochlorothiazide', 'Isotretinoine', 'Tetracycline', 'Hemi-cycline D'] # 化合物の描画 Draw.MolsToGridImage(mols, molsPerRow=2, subImgSize=(450,150), legends=mol_names)
化合物記述子の計算
化合物の比較を行うために1Dと2Dの化合物記述子を抽出、生成します。2D記述子については、あとで化合物の類似度の計算に使用するので、異なるタイプのフィンガープリントを生成します。
1D 化合物記述子:分子量
例示構造の分子量を計算します。
# 化合物の分子量を計算 mol_weights = [Descriptors.MolWt(mol) for mol in mols]
視覚的に比較するために、類似の分子量の化合物構造を描画します。分子量は化合物の類似度にとって有用な記述子となるでしょうか?
# 結果を格納するデータフレームの生成 sim_mw_df = pd.DataFrame({'smiles': smiles, 'name': mol_names, 'mw': mol_weights, "Mol": mols}) # 分子量でソート sim_mw_df.sort_values(['mw'], ascending=False, inplace=True) sim_mw_df[["smiles", "name", "mw"]]
smiles | name | mw | |
---|---|---|---|
0 | CC1C2C(C3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4... | Doxycycline | 444.440 |
6 | CC1(C2CC3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4... | Tetracycline | 444.440 |
3 | CCCCCCCCCCCC(=O)OCCOC(=O)CCCCCCCCCCC | Glycol dilaurate | 426.682 |
1 | CC1(C(N2C(S1)C(C2=O)NC(=O)C(C3=CC=C(C=C3)O)N)C... | Amoxicilline | 365.411 |
2 | C1=COC(=C1)CNC2=CC(=C(C=C2C(=O)O)S(=O)(=O)N)Cl | Furosemide | 330.749 |
5 | CC1=C(C(CCC1)(C)C)C=CC(=CC=CC(=CC(=O)O)C)C | Isotretinoine | 300.442 |
4 | C1NC2=CC(=C(C=C2S(=O)(=O)N1)S(=O)(=O)N)Cl | Hydrochlorothiazide | 297.745 |
7 | CC1C(CC(=O)C2=C1C=CC=C2O)C(=O)O | Hemi-cycline D | 220.224 |
# 分子量とともに化合物を描画 Draw.MolsToGridImage(sim_mw_df["Mol"], legends=[i+': '+str(round(j, 2))+" Da" for i,j in zip(sim_mw_df["name"], sim_mw_df["mw"])], molsPerRow=2, subImgSize=(450, 150))
見てわかるように類似の分子量を持つ化合物は類似の構造をもつことがあります(例 Doxycycline/Tetracycline)。一方で、類似の数の原子を持ちながらも全く異なる原子の配置となっているものもあります(例 Doxycycline/Glycol dilaurate あるいはHydrochlorothiazide/Isotretinoine)。
次に、より詳細な分子の特徴を説明するために、2D化合物記述子を見て見ましょう。
2D 化合物記述子:MACCSフィンガープリント
MACCSフィンガープリントはRDKitを使って簡単に生成することができます。明示的なビットベクトル(explicit bitvector)は我々人間が読めるものではないので、さらにビット文字配列(bitstring)へと変換します。
# MACCSフィンガープリントの生成 maccs_fp1 = MACCSkeys.GenMACCSKeys(mols[0]) # Doxycycline maccs_fp2 = MACCSkeys.GenMACCSKeys(mols[1]) # Amoxicilline maccs_fp1 # <rdkit.DataStructs.cDataStructs.ExplicitBitVect at 0x11baf0ee0>
# フィンガープリントをビット文字配列としてプリント maccs_fp1.ToBitString() # '00000000000000000000000000100000000000000000000000100110000000000010000010101000000011100100110101010100010000101100010000100001000101001001111111101111101011111111110'
# 全化合物のMACCS fingerprintsを生成 maccs_fp_list = [] for i in range(len(mols)): maccs_fp_list.append(MACCSkeys.GenMACCSKeys(mols[i]))
2D 化合物記述子:Morganフィンガープリント
RDKitを使ってMorganサーキュラーフィンガープリントも計算します。2つの異なる関数により、Morganフィンガープリントは整数(int)あるいはビット(bit)ベクトルとして計算することができます。
# Morganフィンガープリントを生成(int vector)、デフォルトでは半径2でベクトルの長さは2048 circ_fp1 = rdFingerprintGenerator.GetCountFPs(mols[:1])[0] circ_fp1 # <rdkit.DataStructs.cDataStructs.UIntSparseIntVect at 0x11badf088>
# セットされた値をみてみます: circ_fp1.GetNonzeroElements() # {45: 1, # 118: 1, # 140: 1, # 163: 1, # 276: 1, # 303: 1, # 309: 1, # 314: 2, # 371: 1, # 438: 1, # 525: 1, # 557: 1, # 650: 3, # 673: 1, # 699: 1, # 807: 6, # 824: 1, # 829: 1, # 881: 1, # 1009: 1, # 1019: 5, # 1027: 1, # 1039: 1, # 1057: 3, # 1060: 1, # 1061: 1, # 1070: 1, # 1082: 1, # 1088: 1, # 1119: 1, # 1154: 1, # 1163: 2, # 1171: 1, # 1257: 1, # 1296: 1, # 1309: 1, # 1341: 1, # 1380: 9, # 1389: 1, # 1457: 1, # 1471: 1, # 1487: 1, # 1582: 1, # 1602: 3, # 1607: 1, # 1630: 1, # 1747: 1, # 1750: 2, # 1831: 1, # 1833: 1, # 1857: 1, # 1873: 3, # 1917: 1, # 1932: 1, # 2000: 1, # 2029: 1}
# Morganフィンガープリントを(bit vectorとして)生成、デフォルトでは半径2でフィンガープリントの長さは2048 circ_b_fp1 = rdFingerprintGenerator.GetFPs(mols[:1])[0] circ_b_fp1 # <rdkit.DataStructs.cDataStructs.ExplicitBitVect at 0x11badf0e0>
# フィンガープリントをビット文字列としてプリント circ_b_fp1.ToBitString() # '00000000000010000000000000000000000000000000000000001000000010000000000000000000000000000000000000000000000000000000000000000000000011000000000010001000000000000000000000000000000011001100000000000000000000000000000000000000000000001000000000000000000000001000000011000000000000000000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000000000000000000000000000000000000000000001000000000000000000000000000000011001100100000000000000000000000000010000000000000000000000000000000000000000000000000000000100000001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000000001100000000000000000000000000000000000000000000000000000000000000000000000000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000100000000000000010000000000000000000000000000000000010000000000000000000000000000000000000001110000010000000000000000000000010000000000000000000000010000000000010000000110000000000110000000000000010000000000000000000000000000000000000000000000000000000000000001100000000000000000000000000000000000000000000000000000000000000000000000000111100000000000000000000000000000000100000000000000010000000100000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000000000000000010000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000001000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000000000000100000000000000000000000000000000000000011100000000000000000'
# 全化合物のMorganフィンガープリントを生成
circ_fp_list = rdFingerprintGenerator.GetFPs(mols)
化合物類似度の計算
次では、2つの類似度評価、すなわちTanimotoとDiceを、2つのタイプのフィンガープリント、すなわちMACCSとMorganフィンガープリントに適用します。
例:2つのMACCSフィンガープリントをタニモト類似度で比較
# 2つの化合物のタニモト係数を計算 DataStructs.TanimotoSimilarity(maccs_fp1, maccs_fp2) # 0.5909090909090909
# 同じ化合物のタニモト係数を計算 DataStructs.TanimotoSimilarity(maccs_fp1, maccs_fp1) # 1.0
次に、クエリ化合物を我々の化合物リストと比較したいと思います。
そこで、RDKitの BulkTanimotoSimilarity
関数とBulkDiceSimilarity
関数を使って、タニモトあるいはDice類似度の類似度評価に基づいて、クエリのフィンガープリントと、リストに格納されたフィンガープリントとの類似度を計算します。
類似度を計算したあとで、次の関数を使ってランク付けした化合物を描画したいと思います。:
def draw_ranked_molecules(sim_df_sorted, sorted_column): """ (ソートした)データフレームの分子を描画する関数 """ # ラベルを定義:最初の分子はクエリ(Query)で、次の分子はランク1から始まる rank = ["#"+str(i)+": " for i in range(0, len(sim_df_sorted))] rank[0] = "Query: " # Doxycyclineと最も類似した化合物(Tanimoto と MACCS フィンガープリント) top_smiles = sim_df_sorted["smiles"].tolist() top_mols = [Chem.MolFromSmiles(i) for i in top_smiles] top_names = [i+j+" ("+str(round(k, 2))+")" for i, j, k in zip(rank, sim_df_sorted["name"].tolist(), sim_df_sorted[sorted_column])] return Draw.MolsToGridImage(top_mols, legends=top_names, molsPerRow=2, subImgSize=(450, 150))
次に、タニモト/Dice類似度評価に基づいて、MACCS/Morganフィンガープリントの比較の全ての組み合わせを調べます。そこで、結果を要約するデータフレームを作成します。
# 結果を格納するデータフレームの生成 sim_df = pd.DataFrame({'smiles': smiles, 'name': mol_names})
MACCSフィンガープリント:タニモト類似度
# 類似度評価スコアをデータフレームに追加 sim_df['tanimoto_MACCS'] = DataStructs.BulkTanimotoSimilarity(maccs_fp1,maccs_fp_list)
# MACCSフィンガープリントのタニモト類似度で並べ替えたデータフレーム sim_df_sorted_t_ma = sim_df.copy() sim_df_sorted_t_ma.sort_values(['tanimoto_MACCS'], ascending=False, inplace=True) sim_df_sorted_t_ma
smiles | name | tanimoto_MACCS | |
---|---|---|---|
0 | CC1C2C(C3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4... | Doxycycline | 1.000000 |
6 | CC1(C2CC3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4... | Tetracycline | 0.928571 |
1 | CC1(C(N2C(S1)C(C2=O)NC(=O)C(C3=CC=C(C=C3)O)N)C... | Amoxicilline | 0.590909 |
7 | CC1C(CC(=O)C2=C1C=CC=C2O)C(=O)O | Hemi-cycline D | 0.403509 |
2 | C1=COC(=C1)CNC2=CC(=C(C=C2C(=O)O)S(=O)(=O)N)Cl | Furosemide | 0.321839 |
4 | C1NC2=CC(=C(C=C2S(=O)(=O)N1)S(=O)(=O)N)Cl | Hydrochlorothiazide | 0.306818 |
5 | CC1=C(C(CCC1)(C)C)C=CC(=CC=CC(=CC(=O)O)C)C | Isotretinoine | 0.288136 |
3 | CCCCCCCCCCCC(=O)OCCOC(=O)CCCCCCCCCCC | Glycol dilaurate | 0.149254 |
# MACCSフィンガープリントのタニモト類似度でランクした分子の描画 draw_ranked_molecules(sim_df_sorted_t_ma, "tanimoto_MACCS")
MACCSフィンガープリントを使用した場合、Tetracyclineは最も類似した分子(スコアが高い)で、ついでAmoxicillineでした。1D 化合物記述子の分子量とは対照的に、線形分子のGlucol dilaurateは類似していない(ランクが最も低い)と認識されました。
MACCSフィンガープリント:Dice類似度
# データフレームへの類似度スコアの追加 sim_df['dice_MACCS'] = DataStructs.BulkDiceSimilarity(maccs_fp1, maccs_fp_list)
# MACCSフィンガープリントのDice類似度でソートしたデータフレーム sim_df_sorted_d_ma = sim_df.copy() sim_df_sorted_d_ma.sort_values(['dice_MACCS'], ascending=False, inplace=True) sim_df_sorted_d_ma
smiles | name | tanimoto_MACCS | dice_MACCS | |
---|---|---|---|---|
0 | CC1C2C(C3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4... | Doxycycline | 1.000000 | 1.000000 |
6 | CC1(C2CC3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4... | Tetracycline | 0.928571 | 0.962963 |
1 | CC1(C(N2C(S1)C(C2=O)NC(=O)C(C3=CC=C(C=C3)O)N)C... | Amoxicilline | 0.590909 | 0.742857 |
7 | CC1C(CC(=O)C2=C1C=CC=C2O)C(=O)O | Hemi-cycline D | 0.403509 | 0.575000 |
2 | C1=COC(=C1)CNC2=CC(=C(C=C2C(=O)O)S(=O)(=O)N)Cl | Furosemide | 0.321839 | 0.486957 |
4 | C1NC2=CC(=C(C=C2S(=O)(=O)N1)S(=O)(=O)N)Cl | Hydrochlorothiazide | 0.306818 | 0.469565 |
5 | CC1=C(C(CCC1)(C)C)C=CC(=CC=CC(=CC(=O)O)C)C | Isotretinoine | 0.288136 | 0.447368 |
3 | CCCCCCCCCCCC(=O)OCCOC(=O)CCCCCCCCCCC | Glycol dilaurate | 0.149254 | 0.259740 |
定義より、タニモトとDice類似度評価は同じランキング結果になりますが、Dice類似度の値の方が大きくなります(タニモトとDiceを求める式はこのトークトリアルの理論編を参照してください)。
Morganフィンガープリント:タニモト類似度
# データフレームへの類似度スコアの追加 sim_df['tanimoto_morgan'] = DataStructs.BulkTanimotoSimilarity(circ_b_fp1, circ_fp_list) sim_df['dice_morgan'] = DataStructs.BulkDiceSimilarity(circ_b_fp1, circ_fp_list)
# Morganフィンガープリントのタニモト類似度で並べ替えたデータフレーム sim_df_sorted_t_mo = sim_df.copy() sim_df_sorted_t_mo.sort_values(['tanimoto_morgan'], ascending=False, inplace=True) sim_df_sorted_t_mo
smiles | name | tanimoto_MACCS | dice_MACCS | tanimoto_morgan | dice_morgan | |
---|---|---|---|---|---|---|
0 | CC1C2C(C3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4... | Doxycycline | 1.000000 | 1.000000 | 1.000000 | 1.000000 |
6 | CC1(C2CC3C(C(=O)C(=C(C3(C(=O)C2=C(C4=C1C=CC=C4... | Tetracycline | 0.928571 | 0.962963 | 0.616279 | 0.762590 |
7 | CC1C(CC(=O)C2=C1C=CC=C2O)C(=O)O | Hemi-cycline D | 0.403509 | 0.575000 | 0.345679 | 0.513761 |
1 | CC1(C(N2C(S1)C(C2=O)NC(=O)C(C3=CC=C(C=C3)O)N)C... | Amoxicilline | 0.590909 | 0.742857 | 0.262136 | 0.415385 |
2 | C1=COC(=C1)CNC2=CC(=C(C=C2C(=O)O)S(=O)(=O)N)Cl | Furosemide | 0.321839 | 0.486957 | 0.141509 | 0.247934 |
5 | CC1=C(C(CCC1)(C)C)C=CC(=CC=CC(=CC(=O)O)C)C | Isotretinoine | 0.288136 | 0.447368 | 0.121495 | 0.216667 |
4 | C1NC2=CC(=C(C=C2S(=O)(=O)N1)S(=O)(=O)N)Cl | Hydrochlorothiazide | 0.306818 | 0.469565 | 0.108911 | 0.196429 |
3 | CCCCCCCCCCCC(=O)OCCOC(=O)CCCCCCCCCCC | Glycol dilaurate | 0.149254 | 0.259740 | 0.084906 | 0.156522 |
# Morganフィンガープリントのタニモト類似度による化合物ランキングの描画 draw_ranked_molecules(sim_df_sorted_t_mo, "tanimoto_morgan")
MACCSとMorganの類似度をタニモト(Morgan) vs タニモト(MACCS)でプロットし比較します。
fig, axes = plt.subplots(figsize=(6,6), nrows=1, ncols=1) sim_df_sorted_t_mo.plot('tanimoto_MACCS','tanimoto_morgan',kind='scatter',ax=axes) plt.plot([0,1],[0,1],'k--') axes.set_xlabel("MACCS") axes.set_ylabel("Morgan") plt.show()
異なるフィンガープリント(ここでは、MACCSフィンガープリントとMorganフィンガープリント)を用いると、異なる類似度の値(ここでは、タニモト係数)となり、ここで示したように、潜在的には化合物類似度のランキングが異なるものとなります。
MorganフィンガープリントはDoxycyclineに対してTetracyclineを(スコアはより低かったでしたが)最も類似した化合物として認識し、Glycol dilaurateを最も似ていない化合物として認識しました。一方で、2番目にランク付されたのはHemi-cycline Dでした。この化合物はサイクリン系化合物の部分構造で、Morganフィンガープリントのアルゴリズムが原子の環境に基づくものであることがその理由であるかもしれません(一方で、MACCSフィンガープリントは特定の特徴の出現頻度を求めるものとなっています)。
類似度検索を使ったバーチャルスクリーニング
フィンガープリントと類似度の計算方法を学んだので、この知識を化合物セット全体からのクエリ化合物の類似度検索に応用することができます。
既知のEGFR阻害剤ゲフィチニブ(Gefitinib)をクエリとして使用し、EGFRに対して試験済みの化合物データセットの中から類似した化合物を検索します。このデータセットは トークトリアル1でChEMBLから抽出し、トークトリアル2でリピンスキーのルールオブファイブによりフィルタリングをおこなったものです。
クエリ化合物をデータセットの全化合物と比較
トークトリアル2で取得したChEMBLデータベースから取り出したEGFRに対して評価済みのフィルタリングされた化合物を含むcsvファイルから化合物を読み込みます。1つのクエリ化合物(ここではゲフィチニブ)を使って、類似の化合物をデータセットの中から探し出します。
# SMILES形式の化合物を含むcsvファイルからデータを読み込む filtered_df = pd.read_csv('../data/T2/EGFR_compounds_lipinski.csv', delimiter=';', usecols=['molecule_chembl_id', 'smiles', 'pIC50']) filtered_df.head()
molecule_chembl_id | smiles | pIC50 | |
---|---|---|---|
0 | CHEMBL63786 | Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 | 11.522879 |
1 | CHEMBL53711 | CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 | 11.221849 |
2 | CHEMBL35820 | CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC | 11.221849 |
3 | CHEMBL53753 | CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 | 11.096910 |
4 | CHEMBL66031 | Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 | 11.096910 |
# クエリ化合物のSMILESからMolオブジェクトを生成 query = Chem.MolFromSmiles('COC1=C(OCCCN2CCOCC2)C=C2C(NC3=CC(Cl)=C(F)C=C3)=NC=NC2=C1'); # Gefitinib, Iressa query
# クエリ化合物のMACCSフィンガープリントとMorganフィンガープリントを生成 maccs_fp_query = MACCSkeys.GenMACCSKeys(query) circ_fp_query = rdFingerprintGenerator.GetCountFPs([query])[0]
# ファイルの全化合物のMACCSフィンガープリントとMorganフィンガープリントを生成 ms = [Chem.MolFromSmiles(i) for i in filtered_df.smiles] circ_fp_list = rdFingerprintGenerator.GetCountFPs(ms) maccs_fp_list = [MACCSkeys.GenMACCSKeys(m) for m in ms]
# クエリ化合物(Gefitinib)とファイルの全化合物のタニモト類似性を計算(MACCS、Morgan)
tanimoto_maccs = DataStructs.BulkTanimotoSimilarity(maccs_fp_query,maccs_fp_list)
tanimoto_circ = DataStructs.BulkTanimotoSimilarity(circ_fp_query,circ_fp_list)
# クエリ化合物(Gefitinib)とファイルの全化合物のDice類似性を計算(MACCS、Morgan)
dice_maccs = DataStructs.BulkDiceSimilarity(maccs_fp_query,maccs_fp_list)
dice_circ = DataStructs.BulkDiceSimilarity(circ_fp_query,circ_fp_list)
# ChEMBL IDとSMILES、Gefitinibに対する化合物のタニモト類似性のテーブルを作成 similarity_df = pd.DataFrame({'ChEMBL_ID':filtered_df.molecule_chembl_id, 'bioactivity':filtered_df.pIC50, 'tanimoto_MACCS': tanimoto_maccs, 'tanimoto_morgan': tanimoto_circ, 'dice_MACCS': dice_maccs, 'dice_morgan': dice_circ, 'smiles': filtered_df.smiles,})
# データフレームを表示
similarity_df.head()
ChEMBL_ID | bioactivity | tanimoto_MACCS | tanimoto_morgan | dice_MACCS | dice_morgan | smiles | |
---|---|---|---|---|---|---|---|
0 | CHEMBL63786 | 11.522879 | 0.409836 | 0.324786 | 0.581395 | 0.490323 | Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 |
1 | CHEMBL53711 | 11.221849 | 0.484375 | 0.327434 | 0.652632 | 0.493333 | CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 |
2 | CHEMBL35820 | 11.221849 | 0.666667 | 0.445455 | 0.800000 | 0.616352 | CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC |
3 | CHEMBL53753 | 11.096910 | 0.428571 | 0.333333 | 0.600000 | 0.500000 | CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 |
4 | CHEMBL66031 | 11.096910 | 0.384615 | 0.345133 | 0.555556 | 0.513158 | Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 |
類似性評価の値の分布
理論編で述べたように、同じフィンガープリント(例 MACCSフィンガープリント)について比較すれば、タニモト類似度の値はDIce類似度の値よりも小さくなります。また、2つの異なるフィンガープリント(例 MACCSフィンガープリントとMorganフィンガープリント)を比較すると、類似性評価の値(例 タニモト類似度)は変化します。
ヒストグラムをプロットすることで分布を見ることができます。
# MACCSフィンガープリントのタニモト類似度の分布をプロット %matplotlib inline fig, axes = plt.subplots(figsize=(10,6), nrows=2, ncols=2) similarity_df.hist(["tanimoto_MACCS"], ax=axes[0,0]) similarity_df.hist(["tanimoto_morgan"], ax=axes[0,1]) similarity_df.hist(["dice_MACCS"], ax=axes[1,0]) similarity_df.hist(["dice_morgan"], ax=axes[1,1]) axes[1,0].set_xlabel("similarity value") axes[1,0].set_ylabel("# molecules") plt.show()
ここでも類似度を比較します。今回は直接、2つのフィンガープリントに関するタニモト類似度とDice類似度を比較します。
fig, axes = plt.subplots(figsize=(12,6), nrows=1, ncols=2) similarity_df.plot('tanimoto_MACCS','dice_MACCS',kind='scatter',ax=axes[0]) axes[0].plot([0,1],[0,1],'k--') axes[0].set_xlabel("Tanimoto(MACCS)") axes[0].set_ylabel("Dice(MACCS)") similarity_df.plot('tanimoto_morgan','dice_morgan',kind='scatter',ax=axes[1]) axes[1].plot([0,1],[0,1],'k--') axes[1].set_xlabel("Tanimoto(Morgan)") axes[1].set_ylabel("Dice(Morgan)") plt.show()
類似度分布は類似度を解釈するのに重要です(例 MACCSフィンガープリントとMorganフィンガープリント、タニモト類似度とDice類似度について値0.6は異なる評価を与えられる必要があります)
次では、Morganフィンガープリントに基づき、タニモト類似度で最もよく似た化合物を描画します。
最も類似の化合物を描画
私たちの作成したランキングにおいて最も類似した化合物との比較として、ゲフィチニブ(Gefitinib)の構造を視覚的に調べます。生理活性の情報(トークトリアル1でChEMBLから抽出したpIC50)も含めます。
# tanimoto_morganでソートしたデータフレーム similarity_df.sort_values(['tanimoto_morgan'], ascending=False, inplace=True) similarity_df.head()
ChEMBL_ID | bioactivity | tanimoto_MACCS | tanimoto_morgan | dice_MACCS | dice_morgan | smiles | |
---|---|---|---|---|---|---|---|
2563 | CHEMBL939 | 6.288193 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCOCC1 |
877 | CHEMBL14699 | 8.000000 | 1.000000 | 0.923913 | 1.000000 | 0.960452 | COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCN1CCOCC1 |
1693 | CHEMBL299672 | 7.148742 | 0.919355 | 0.843750 | 0.957983 | 0.915254 | COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCCC1 |
1867 | CHEMBL24137 | 7.000000 | 1.000000 | 0.836735 | 1.000000 | 0.911111 | COc1cc2c(Nc3ccc(Cl)cc3F)ncnc2cc1OCCCN1CCOCC1 |
1544 | CHEMBL384699 | 7.292430 | 1.000000 | 0.836735 | 1.000000 | 0.911111 | COc1cc2ncnc(Nc3cc(Cl)ccc3F)c2cc1OCCCN1CCOCC1 |
# データフレームにSMILES文字列の構造表現(ROMol - RDKit オブジェクト Mol)を追加 PandasTools.AddMoleculeColumnToFrame(similarity_df, 'smiles')
# クエリ構造とトップランクの化合物群(+生理活性)の描画 sim_mols = [Chem.MolFromSmiles(i) for i in similarity_df.smiles][:11] legend = ['#' + str(a) + ' ' + b + ' ('+str(round(c,2))+')' for a, b, c in zip(range(1,len(sim_mols)+1), similarity_df.ChEMBL_ID, similarity_df.bioactivity)] Chem.Draw.MolsToGridImage(mols = [query] + sim_mols[:11], legends = (['Gefitinib'] + legend), molsPerRow = 4)
データセットにおいてゲフィチニブと比較してトップにランクした化合物群は、最初は我々のデータセットに含まれるゲフィチニブのエントリー(rank1 と rank2)で、続いてゲフィチニブの変換体(例 benzole置換基パターンが異なるもの)です。
注:ChEMBLにはゲフィチニブ(よく研究された化合物なので)の完全な構造活性相関分析がふくまれていて、したがって私たちが取得したデータセットにゲフィチニブ様化合物が多く含まれていることは驚くべきことではありません。
それでは、類似度検索がどの程度、データセット上の活性化合物と不活性化合物を区別することができるか、その性能をチェックしたいと思います。そこで、トークトリアル1 でChEMBLから取得した化合物の(EGFRに対する)生理活性の値を使用します。
エンリッチメントプロットの生成
バーチャルスクリーニングの妥当性を評価し、見つかった活性化合物の比率を見るためにエンリッチメントプロットを作成します。
エンリッチメントプロットが示すのは; * データセット全体のうち、トップにランクした化合物の比率(x-axis) vs. * データセット全体のうち活性化合物(y-axis)の比率
MACCSフィンガープリントとMorganフィンガープリントのタニモト類似度を比較します。
化合物を活性化合物あるいは不活性化合物のいずれとして取り扱うかを決めるために、一般に使用されるpIC50のカットオフ値6.3を適用します。文献中にはpIC50カットオフ値として5〜7にわたる範囲でいくつか提案がなされていて、データポイントをとらない排除範囲を定義しているものもありますが、私たちはこのカットオフ(6.3)は合理的と考えています。 同じカットオフをトークトリアル10の機械学習でも用います。
# 活性化合物と不活性化合物を区別するpIC50 カットオフ値 threshold = 6.3
similarity_df.head()
ChEMBL_ID | bioactivity | tanimoto_MACCS | tanimoto_morgan | dice_MACCS | dice_morgan | smiles | |
---|---|---|---|---|---|---|---|
2563 | CHEMBL939 | 6.288193 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCOCC1 |
877 | CHEMBL14699 | 8.000000 | 1.000000 | 0.923913 | 1.000000 | 0.960452 | COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCN1CCOCC1 |
1693 | CHEMBL299672 | 7.148742 | 0.919355 | 0.843750 | 0.957983 | 0.915254 | COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCCC1 |
1867 | CHEMBL24137 | 7.000000 | 1.000000 | 0.836735 | 1.000000 | 0.911111 | COc1cc2c(Nc3ccc(Cl)cc3F)ncnc2cc1OCCCN1CCOCC1 |
1544 | CHEMBL384699 | 7.292430 | 1.000000 | 0.836735 | 1.000000 | 0.911111 | COc1cc2ncnc(Nc3cc(Cl)ccc3F)c2cc1OCCCN1CCOCC1 |
def get_enrichment_data(similarity_df, similarity_measure, threshold): """ エンリッチメントプロットのxとyの値を計算する関数: x - データセットで何%にランクしているか y - 何%の本当に活性な化合物が見つかったか """ # データセットの化合物の数を取得 mols_all = len(similarity_df) # データセットの活性化合物の数を取得 actives_all = sum(similarity_df.bioactivity >= threshold) # データセット全体を処理している間、活性化合物のカウンターを保持するリストを初期化 actives_counter_list = [] # 活性化合物のカウンターを初期化 actives_counter = 0 # 注: エンリッチメントプロットのためデータをランク付けしなければなりません。 # 選択した類似度評価によって化合物を並べ替えます。 similarity_df.sort_values([similarity_measure], ascending=False, inplace=True) # ランク付けされたデータセットを一つずつ処理し、(生理活性をチェックすることで)各化合物が活性化合物どうか確認します for value in similarity_df.bioactivity: if value >= threshold: actives_counter += 1 actives_counter_list.append(actives_counter) # 化合物の数をデータセットのランク何%になるかに変換 mols_perc_list = [i/mols_all for i in list(range(1, mols_all+1))] # 活性化合物の数を本当の活性化合物の何%が見つかったかに変換 actives_perc_list = [i/actives_all for i in actives_counter_list] # xとyの値とラベルをもつデータフレームを生成 enrich_df = pd.DataFrame({'% ranked dataset':mols_perc_list, '% true actives identified':actives_perc_list, 'similarity_measure': similarity_measure}) return enrich_df
# プロットする類似度評価を定義 sim_measures = ['tanimoto_MACCS', 'tanimoto_morgan'] # 全類似度評価についてエンリッチメントプロットのデータを持つデータフレームのリストを作成 enrich_data = [get_enrichment_data(similarity_df, i, threshold) for i in sim_measures]
# プロットのためのデータセットを準備: # 類似度評価毎のデータフレームを一つのデータフレームに連結 # …異なる類似度評価は「similarity_measure」列によって区別可能です enrich_df = pd.concat(enrich_data)
fig, ax = plt.subplots(figsize=(6, 6)) fontsize = 20 for key, grp in enrich_df.groupby(['similarity_measure']): ax = grp.plot(ax = ax, x = '% ranked dataset', y = '% true actives identified', label=key, alpha=0.5, linewidth=4) ax.set_ylabel('% True actives identified', size=fontsize) ax.set_xlabel('% Ranked dataset', size=fontsize) # データセットの活性化合物比率 ratio = sum(similarity_df.bioactivity >= threshold) / len(similarity_df) # 理想的な場合のカーブをプロット ax.plot([0,ratio,1], [0,1,1], label="Optimal curve", color="black", linestyle="--") # ランダムな場合のカーブをプロット ax.plot([0,1], [0,1], label="Random curve", color="grey", linestyle="--") plt.tick_params(labelsize=16) plt.legend(labels=['MACCS', 'Morgan', "Optimal", "Random"], loc=(.5, 0.08), fontsize=fontsize, labelspacing=0.3) # プロットを保存ーテキストボックスを含めるためにbbox_inchesを使用: # https://stackoverflow.com/questions/44642082/text-or-legend-cut-from-matplotlib-figure-on-savefig?rq=1 plt.savefig("../data/T4/enrichment_plot.png", dpi=300, bbox_inches="tight", transparent=True) plt.show()
エンリッチメントプロットによるとMACCSフィンガープリントよりもMorganフィンガープリント基づく比較の方が少し良いパフォーマンスを示しています。
# ランク付されたデータセットのx%についてEFを取得 def print_data_ef(perc_ranked_dataset, enrich_df): data_ef = enrich_df[enrich_df['% ranked dataset'] <= perc_ranked_dataset].tail(1) data_ef = round(float(data_ef['% true actives identified']), 1) print("Experimental EF for ", perc_ranked_dataset, "% of ranked dataset: ", data_ef, "%", sep="") # ランク付されたデータセットのx%についてランダムEFを取得 def print_random_ef(perc_ranked_dataset): random_ef = round(float(perc_ranked_dataset), 1) print("Random EF for ", perc_ranked_dataset, "% of ranked dataset: ", random_ef, "%", sep="") # ランク付されたデータセットのx%について理想的な場合のEFを取得 def print_optimal_ef(perc_ranked_dataset, similarity_df, threshold): ratio = sum(similarity_df.bioactivity >= threshold) / len(similarity_df) * 100 if perc_ranked_dataset <= ratio: optimal_ef = round(100/ratio * perc_ranked_dataset, 1) else: optimal_ef = round(float(100), 2) print("Optimal EF for ", perc_ranked_dataset, "% of ranked dataset: ", optimal_ef, "%", sep="")
# パーセンテージを選択 perc_ranked_list = 5 # EFデータを取得 print_data_ef(perc_ranked_list, enrich_df) print_random_ef(perc_ranked_list) print_optimal_ef(perc_ranked_list, similarity_df, threshold) # Experimental EF for 5% of ranked dataset: 1.0% # Random EF for 5% of ranked dataset: 5.0% # Optimal EF for 5% of ranked dataset: 8.9%
訳注(04/2020)
オリジナルの実践編はここまでですが、このEFの結果は少しおかしい気がします。エンリッチメントプロットを見ると、エンリッチメントファクターは「optimal > Experimental > Random」となると思われます。RandomよりもExperimentalが低いということは、むしろ不活性化合物へと選択のバイアスがかかっていることになってしまいます。
どこかおかしいところがないか?順番に見ていきます。
まずEFの算出に使われているDataFrame enrich_dfは2つの類似度評価基準のデータを繋げたものなので、それぞれ別々にしてみます。
# tanimoto_MACCS enrich_df_taMA = enrich_df[enrich_df['similarity_measure'] == 'tanimoto_MACCS'] # tanimoto_morgan enrich_df_tamo = enrich_df[enrich_df['similarity_measure'] == 'tanimoto_morgan'] print("Size of enrich_df: ", len(enrich_df)) print("Size of tanimoto_MACCS DataFrame: ", len(enrich_df_taMA)) print("Size of tanimoto_morgan DataFrame: ", len(enrich_df_tamo)) # Size of enrich_df: 9046 # Size of tanimoto_MACCS DataFrame: 4523 # Size of tanimoto_morgan DataFrame: 4523
DataFrameの5% ranked Datasetに相当する箇所を見てみます。
# 5% に相当する数 index_5perc = round(len(enrich_df_taMA)*0.05) # DataFrameのindexは0から始まるので-1した行を表示 enrich_df_taMA[index_5perc-1:index_5perc]
% ranked dataset | % true actives identified | similarity_measure | |
---|---|---|---|
225 | 0.049967 | 0.07319 | tanimoto_MACCS |
見やすさのためsliceでDataFrameとして取り出しています。
ランク上位5%(0.049)に相当する数のなかに、実際の活性評価でactiveだったものは7.3%(% true actives identified, 0.07319)となっています。この値は先のRandom、Optimalと比較して妥当な値に思います。
DataFrameのデータ自体には問題なさそうなので、値の取り出し方(関数print_data_ef
)に問題があったのでしょう。
関数の中身を順番に実行してみます。
# 5%に設定 perc_ranked_dataset = 5 # 5%内のDataFrameを取り出し、その一番最後の行(tail)を取り出す enrich_df[enrich_df['% ranked dataset'] <= perc_ranked_dataset].tail(1)
% ranked dataset | % true actives identified | similarity_measure | |
---|---|---|---|
4522 | 1.0 | 1.0 | tanimoto_morgan |
取り出されたのは「index:4522」で4523番目の行です。この% true actives identifed
列がEFとして取り出されていた値(1.0)です。
閾値5以下で取り出されたのはsimilarity_measure
:tanimoto_morganの全化合物でした。単純にDataFrameのデータは%に換算していないのに、取り出す際に%換算の値を使ってしまったのが原因のようです。
それではそれぞれの類似度について正しい値を確認してみます。
# 関数の再定義 def print_data_ef2(perc_ranked_dataset, enrich_df): perc_ranked_dataset_100 = perc_ranked_dataset / 100 data_ef = enrich_df[enrich_df['% ranked dataset'] <= perc_ranked_dataset_100].tail(1) data_ef = round(float(data_ef['% true actives identified'] * 100), 1) print("Experimental EF for ", perc_ranked_dataset, "% of ranked dataset: ", data_ef, "%", sep="")
# MACCS keyの場合 # パーセンテージを選択 perc_ranked_list = 5 # EFデータを取得 print_data_ef2(perc_ranked_list, enrich_df_taMA) print_random_ef(perc_ranked_list) print_optimal_ef(perc_ranked_list, similarity_df, threshold) # Experimental EF for 5% of ranked dataset: 7.3% # Random EF for 5% of ranked dataset: 5.0% # Optimal EF for 5% of ranked dataset: 8.9%
# Morganフィンガープリントの場合 # パーセンテージを選択 perc_ranked_list = 5 # EFデータを取得 print_data_ef2(perc_ranked_list, enrich_df_tamo) print_random_ef(perc_ranked_list) print_optimal_ef(perc_ranked_list, similarity_df, threshold) # Experimental EF for 5% of ranked dataset: 7.9% # Random EF for 5% of ranked dataset: 5.0% # Optimal EF for 5% of ranked dataset: 8.9%
いずれも「optimal > Experimental > Random」となっており、Morganの方がMACCSよりも若干良い値となっています。無事エンリッチメントプロットと比較してもおかしくない値が得られました。
訳注ここまで
ディスカッション
ここではタニモト類似度を使ってバーチャルスクリーニングを実施しました。もちろん、Dice類似度や他の類似度評価を使っても行うことができます。
化合物フィンガープリントを使用した類似度検索の欠点は、化合物類似度に基づくものなので新規な構造を生み出さないことです。化合物類似度を扱う上でのその他の課題としては、いわゆるアクティビティクリフ(activity clliff)があります。分子の官能基におけるちょっとした変化が生理活性の大きな変化を起こすことがあります。
クイズ
- アクティビティクリフを回避するにはどこから始めれば良いでしょうか?
- MACCSフィンガープリントとMorganフィンガープリントを互いに比較した場合の利点と欠点は何でしょう?
- 使用したフィンガープリントによっておこる、類似度データフレームにおける順序の違いをどう説明できるでしょうか?
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。
フィンガープリントと分子の類似性について:
* 化学の新しいカタチ:RDKitでフィンガープリントを使った分子類似性の判定
また、同書の7章: グラフ構造を利用した類似性の評価にはActivity Cliffの説明と事例も取り上げられています。全体を通して非常に勉強になるので未読の方は是非一読をお勧めします。
エンリッチメントファクターとその他の指標については以下の記事が参考になりました。
tonetsの日記: バーチャルスクリーニングの指標のテスト その2
難しすぎて私では理解できませんでしたがECFPを実装されている方もいらっしゃいました。
創薬・材料探索のための機械学習: ECFPとNeural-fingerprintの比較
私は既存のメソッドの使い方ですらよくわからん、誰か日本語で解説して!という感じなので実際にいろいろ作ってしまう方々は本当にすごいなあと思います。アホですみません。
勝手にいろいろと引用させていただいておりますので問題があれば削除します。