magattacaのブログ

日付以外誤報

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

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

以下、訳。

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

トークトリアル 10

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

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

Angelika Szengel, Marvis Sydow and Dominique Sydow

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

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

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

学習の目標

理論

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

実践

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

レファレンス

結合サイト比較:

イマチニブ(Imatinib):

理論

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

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

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

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

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

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

f:id:magattaca:20200506005325p:plain

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

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

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

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

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

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

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

実践

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

import os
import pprint
import pickle
import glob
import time

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

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

from pypdb import *
from biopandas.pdb import PandasPdb

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

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

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

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

f:id:magattaca:20200506005407p:plain

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

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

sti_mol = Chem.AddHs(sti)
sti_mol

f:id:magattaca:20200506005426p:plain

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

f:id:magattaca:20200506005445p:plain

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

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

os.popen('pymol -R')

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

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

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

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

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

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

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

objPMV.GetPNG(h=300)

f:id:magattaca:20200506005532p:plain

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

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

PDBを検索

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

PDB構造の可視化

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

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

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

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

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

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

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

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

png

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

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

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

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

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

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

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

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

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

objPMV.GetPNG(h=500)

f:id:magattaca:20200506005836p:plain

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

objPMV.HideObject("3FW1")

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

objPMV.GetPNG(h=500)

f:id:magattaca:20200506005855p:plain

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

f:id:magattaca:20200506010049p:plain

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

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

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

f:id:magattaca:20200506010140p:plain

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

f:id:magattaca:20200506010215p:plain

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

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

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

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

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

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

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

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

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

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

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

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

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

f:id:magattaca:20200506010314p:plain

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

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

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

ディスカッション

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

クイズ

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

訳以上

追記

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

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

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

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

TeachOpenCADD トピック9 〜リガンドベースファーマコフォア〜

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

以下、訳。

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

トークトリアル 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)を抽出
  • 全リガンドのファーマコフォアフィーチャーを表示
    • 水素結合ドナー
    • 水素結合アクセプター
    • 疎水性コンタクト
  • フィーチャータイプ毎にフィーチャーの座標を収集
  • 組み合わせファーマコフォアの作成
  • クラスターの表示
    • 水素結合ドナー
    • 水素結合アクセプター
    • 疎水性コンタクト
  • 組み合わせファーマコフォアの表示

レファレンス

理論

ファーマコフォア

コンピューター支援医薬品デザイン(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)をお勧めします。

f:id:magattaca:20200505003716p:plain:w400

Figure 1: タンパク質-リガンド相互作用を表す構造ベースのファーマコフォア(Dominique Sydowによる図)

ファーマコフォアを用いたバーチャルスクリーニング

先にトークトリアル 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平均法クラスタリングを使います:

  1. k個の異なる重心(centroid)を選択し、データセットの各ポイントを最も近い重心に割り当てる
  2. 現在のクラスターに基づいて新しい重心を計算し、データセットの各ポイントが新しい重心のうち最も近いものに割り当てられる
  3. 重心が安定化するまでこの過程を繰り返す

(K means wikipedia)

実践

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])

f:id:magattaca:20200505003924p:plain

ここで問題に行き当たりました: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])

f:id:magattaca:20200505004018p:plain

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])

f:id:magattaca:20200505004038p:plain

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()

f:id:magattaca:20200505004232p:plain

ファーマコフォアフィーチャーの抽出

上で述べたように、このトークトリアルの目的はリガンドセットからリガンドベースのファーマコフォアの組み合わせを作成することです。まず最初に、リガンド毎にファーマコフォアフィーチャーを取り出す必要があります。そこで、フィーチャーファクトリーを(デフォルトのフィーチャーの定義で)読み込みます。

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標準ライブラリcollectionsCounterクラスにより、与えられたリストの要素の出現頻度について辞書型オブジェクトが生成されています。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()

f:id:magattaca:20200505004708p:plain

水素結合アクセプター

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()

f:id:magattaca:20200505004856p:plain

疎水性コンタクト

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()

f:id:magattaca:20200505004938p:plain

フィーチャータイプ毎にフィーチャーの座標を集める

(フィーチャータイプ毎に)フィーチャーをクラスタリングしたいので、(フィーチャータイプ毎に)フィーチャーのすべての座標を集めます。

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は、collectionsCounterクラスターのラベルを集計(要素の出現頻度の辞書を作成)した後に並べ替えています。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()

f:id:magattaca:20200505005317p:plain

水素結合アクセプター

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()

f:id:magattaca:20200505005357p:plain

疎水性コンタクト

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()

f:id:magattaca:20200505005435p:plain

組み合わせファーマコフォアの表示

この最後のステップでは、クラスタリングしたファーマコフォアのフィーチャー(即ち、水素結合ドナー、アクセプターと疎水性コンタクト)を組み合わせて、一つの組み合わせファーマコフォアを作成します。この組み合わせファーマコフォアは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()

f:id:magattaca:20200505005507p:plain

ディスカッション

このトークトリアルでは、EGFRに結合することが知られているリガンドのセットを、予め重ね合わせたものを用いて、組み合わせファーマコフォアモデルを作成しました。これで、既知のリガンドに観測された立体的、物理化学的特徴を示し、EGFR結合サイトに結合するかもしれないような新奇な低分子を見つけることを目的として、低分子化合物の巨大なライブラリのバーチャルスクリーニングを行うために、この組み合わせファーマコフォアモデルを使うことができます。

スクリーニングの前に、通常、ファーアマコフォアモデルはさらに最適化されます。例えば、スクリーニングに用いるフィーチャーの特徴の数を減らすために、生物学的な知識(ある相互作用は重要である一方、他の相互作用はそうではないという報告されているかもしれません)、あるいは化学的専門知識に基づいて、いくつかのフィーチャーは除外されるかもしれません。

このトークトリアルではバーチャルスクリーニングは扱いませんが、RDKitを使ってファーマコフォアモデリングとバーチャルスクリーニングをデモンストレーションしている、Nik Stieflによる素晴らしいチュートリアルに言及しておきます。(RDKit UGM 2016 on GitHub).

ファーマコフォアフィーチャーのクラスタリングを行うためにk平均法クラスタリングを用いました。このクラスタリング手法は、使用者が予めクラスターの数を定義する必要があるという欠点があります。通常、クラスターの数は、クラスタリングの実施前(あるいはクラスターの精度を高めていく過程)に、ポイントの分布を目視で調査して決めるので、ファーマコフォアを自動的に作成するには妨げとなります。この解決策として、密度に基づくクラスタリング手法(とk平均法クラスタリングの組み合わせも)が挙げられます。

クイズ

  1. ファーマコフォアフィーチャーとファーマコフォアを説明してください。
  2. 構造ベースのファーマコフォアモデリングとリガンドベースのファーマコフォアモデリングの違いを説明してください。
  3. 組み合わせファーマコフォアを導く方法を説明してください。

訳以上

追記

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

今回のトピックに関連して日本語で読める記事として株式会社理論創薬研究所 吉森篤史氏による以下の記事をお勧めします。 RDKitを利用したファーマコフォアの定義とその使用方法について解説がされていてとてもわかりやすいです。

欲を言えば作ったファーマコフォアモデルによるバーチャルスクリーニングの実施まで行なってほしかったですが、計算リソースとデータセットの準備だけでも大変そうですね。

勝手にいろいろと引用させていただいておりますので問題があれば削除します。

TeachOpenCADD トピック8 〜タンパク質データの取得〜

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

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

以下、訳。

トークトリアル 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で重ね合わせ(アラインメント)、次のトークトリアルで使用するためにリガンドを抜き出し保存します。

学習の目標

理論

  • Protein Data Bank (PDB)
  • Pythonパッケージ PyPDB

実践

  • 検索するクエリタンパク質の選択
  • クエリタンパク質のPDBエントリーに関する統計値を取得
  • クエリタンパク質の全PDB IDを取得
  • PDBエントリーのメタ情報を取得
  • PDBエントリーのメタ情報をフィルタリングし並べ替え
  • 上位の構造からリガンドのメタ情報を取得
  • 上位リガンド化合物を描画
  • タンパク質ーリガンド IDのペアを作成
  • PDB構造ファイルの取得
  • PDB構造の重ね合わせ

レファレンス

理論

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 で見ることができます。

f:id:magattaca:20200503195241p:plain:w400

Figure 1: 上皮増殖因子(EGFR)の例としてPDB ID 3UG5のタンパク質構造(灰色)と相互作用するリガンド(緑色)を描画(Dominique Sydowによる図)。

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()

f:id:magattaca:20200503195456p:plain

クエリタンパク質の全PDB IDを取得

次に、pypdb関数のmake_querydo_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)

f:id:magattaca:20200503195823p:plain

タンパク質ーリガンド 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構造ファイルの取得

では、pypdbget_pdb_file関数を使ってPDBからPDB構造ファイルを取得します。すなわちタンパク質、リガンド(そして可能な場合には補因子や水分子、そしてイオンといった他の原子あるいは分子)の3次元座標を取得します。 利用可能なファイルフォーマットはpdbcifで、タンパク質(とリガンド、補因子、水分子、イオン)の原子の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に渡す引数のrefAlignTargetrefAlignQueryはどちらも( 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を参照)と、リガンドの生理活性をチェックする(即ち、活性があるリガンドだけを含める)ことをお勧めします。

クイズ

  1. Protein Data Bankの含むデータの種類を要約してください。
  2. 構造の分解能が何を表し、このトークトリアルではどうやって、また、なぜ分解能についてフィルタリングしたのかを説明してください。
  3. 構造の重ね合わせ(アラインメントを取ること)とはどういう意味か説明し、このトークとリアルで実施したアラインメントについて議論してください。

訳以上

追記

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

今回、 取得した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の取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!

以下、訳。

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

トークトリアル 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阻害剤となる可能性のある化合物を見つけるための、機械学習に基づくスクリーニングパイプラインの構築と評価

レファレンス


理論

f:id:magattaca:20200502110047p:plain:w200,right

機械学習をうまく適用するために必要となるものは、巨大な化合物データセット、化合物のエンコーディング、データセットの化合物ごとのラベル、そしてモデルを訓練するための機械学習アルゴリズムです。これらによって、新しい化合物の予測を行うことができます。

データの準備

機械学習のために、化合物を特徴のリストに変換する必要があります。化合物の表現方法として、フィンガープリントがよく使われます。

このトークトリアルで使うフィンガープリントはRDKitに実装されているものです(より詳細な情報はG. Landrumのプレゼンテーション で見ることができます)。

  • maccs: MACCS キー 166ビットの構造キー記述子(structural key descriptors)で、各ビットはSMARTSパターンに関連づけられています。
  • ecfp4ecfp6: 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 参照)

教師あり学習

学習アルゴリズムは訓練データの中のパターンを見つけることでルールを作成します。

f:id:magattaca:20200502110557p:plain:w250,right

  • ランダムフォレスト(Random Forest、RF):多数の決定木(decision tree)により平均化された予測を出します。
  • サポートベクトルマシン(Support Vector Machines (SVM)): SVMカーネルトリックと呼ばれる入力を高次元の特徴空間に暗にマッピングする手法を使って、非線形の分類を効率的に実行することができます。オブジェクトの関数としてマージンを最大化するというアイデアに基づいた分類です。

f:id:magattaca:20200502110625p:plain:w150,right

  • 人工ニューラルネットワーク(Artificial neural networks、ANNs):ANNは生物学的な脳の神経をゆるくモデル化した人工神経(artificial neuron)と呼ばれる連結されたユニットあるいはノードの集まりに基づいています。生物学的な脳のシナプスの様に、各結合は一つの人工神経からもう一方へとシグナルを伝達することができます。シグナルを受け取った人工神経は、そのシグナルを処理し、結合する別の人工神経へとシグナルを送ります。(図はWikipediaより)

妥当性検査手法(Validation strategy):k分割交差検証

  • このモデル妥当性評価テクニックは反復的な方法でデータセットを2つのグループに分けます。

  • この検証の目標は過学習(over-fitting)として知られる問題に注意するために、モデルが一度も見たことのないデータを予測する性能をテストし、モデルの汎化性能を評価することです。

性能評価

f:id:magattaca:20200502110705p:plain:w250,right

  • 感度(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)の診断性能を図解するプロット図
    • 感度を特異度に対してプロットします。
  • AUCROC曲線の下の面積(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

f:id:magattaca:20200502111529p:plain

我々のモデルは全ての指標について非常に良い値を示し、予測性能があるように見えます。

サポートベクトル分類器

ここでは動径基底関数カーネル(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

f:id:magattaca:20200502111632p:plain

ニューラルネットワーク分類器

ここで試す最後のアプローチはニューラルネットワークモデルです。各層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

f:id:magattaca:20200502111727p:plain

ディスカッション

  • 我々のデータッセットで最もパフォーマンスがよかったのはどのモデルで、それは何故でしょうか?

    • データセットに対して、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が非常に充実していると思います。 取り上げられた内容をまとめてくださっているページもあるので、そちらを参照いただいて関連する記事にアクセスするのが一番早いかと思います。

ざっと見たところ今回の内容と関連するところとしては

などがありそうです。素敵!

勝手にいろいろと引用させていただいておりますので問題があれば削除します。

TeachOpenCADD トピック6 〜最大共通部分構造〜

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

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

以下、訳。

トークトリアル 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アルゴリズムの詳細な説明

実践

レファレンス


理論

導入

f:id:magattaca:20200426180517p:plain:w300

最大共通部分構造(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): O6rdkit 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)

f:id:magattaca:20200426180706p:plain:w350

  • MCSの部分とはなり得ない結合の除去
  • 各入力構造に原子タイプと結合タイプの情報が存在している必要があります
  • 結合タイプ:最初の原子、結合、2つめの原子のSMATSからなる文字列
  • 入力構造全てに存在していない全ての結合タイプを排除し、それぞれの辺(edge)を削除
  • 結果:全ての原子の情報をもちつつ、辺(edge(結合))の数がより少ないフラグメント化された構造

最大フラグメントのうち最小の構造をクエリとして用いる

f:id:magattaca:20200426180815p:plain:w150

  • ヒューリスティックなアプローチ:
    • 各入力構造の最大のフラグメントを見つける
    • 最大フラグメントの結合の数によって入力構造を昇順に並べ換える
    • 原子の数との紐付けを解き、入力の順を代わりのものとして使う
  • 最大フラグメントのうち最小の構造をクエリ構造とする
  • 他の入力構造由来の構造がターゲットとなります

フラグメントのサブグラフを列挙するために、幅優先探索(breadth-first search (BFS))と優先度付きキュー(priority que)を使う

f:id:magattaca:20200426180845p:plain:w150

  • いわゆる種(シード、seed)を成長させることに基づく列挙(enumeration)
  • シード:現在のサブグラフの原子/結合と、除外するセット(成長には使えない結合)
  • 冗長性を防ぐために:
    • 最初のシードはフラグメントの最初の結合で、フラグメント全体のサイズにまで成長する可能性があります
    • 2番目のシードは2番目の結合で、最初の結合を使ったものからは排除されます
    • 3番目のシードは3番目の結合から始まり、最初と2番目を使用したものからは排除されます
    • ・・・
  • シードはシードにつながっている結合に沿って成長します(除外セットや既にシードに含まれているものは含みません)
  • 各ステップですべての成長の可能性が考慮されます
  • 例、拡張の方向としてN個の可能な結合がある時、2のN-1乗個の新しいシードがキューにつけ加えられる

f:id:magattaca:20200426180914p:plain:w350

  • サブグラフに加える新しい結合がなくなった時点で列挙は終わります(シードがキューから除外されます)
  • 最大のシードが最初に処理されます

他の全ての構造に含まれていないシードを取り除く

  • 各成長の段階で → 新しいシードが他の全ての構造に存在しているかどうかチェックする  
  • 存在していない場合:キューからシードを取り除く

十分な成長のポテンシャルがないシードを取り除く

  • 除外リストと拡張の可能性がある辺(エッジ)から成長ポテンシャルを評価します
  • 現在の最も良いサブグラフよりもポテンシャルが小さい場合 -> シードをキューから取り除きます

これらのアプローチを使うことで、最大共通部分構造に相当する最大のサブグラフを順に辿っていくことことは普通の問題となります。

実践

# 必要なモジュールをインポート
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

f:id:magattaca:20200426185241p:plain

様々なパラメーター入力値で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]

f:id:magattaca:20200426185312p:plain

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)

f:id:magattaca:20200426185337p:plain

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]

f:id:magattaca:20200426185404p:plain

highlightMolecules(mols, mcs2, 12)

f:id:magattaca:20200426185420p:plain

この例を見ればわかるように、化合物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'])

f:id:magattaca:20200426185441p:plain

ここでは選択した閾値とパラメーターに依存して、少し異なる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 Mol
1 CHEMBL53711 nM 0.006 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849 Mol
2 CHEMBL35820 nM 0.006 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849 Mol
# 選択した高活性化合物だけについて計算を実行します
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]

f:id:magattaca:20200426185539p:plain

highlightMolecules(rand_mols, mcs2, 12)

f:id:magattaca:20200426185601p:plain

閾値インタラクティブに変更することもできます。

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の取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

以下、訳。

トークトリアル 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)の構築を行うことができます。そのようなクラスタリングによって、より大きなスクリーニング化合物セットから、さらなる実験を行う対象となる化合物の多様性のあるセットを選ぶことができます。

学習の目標

このトークトリアルでは、以下についてさらに学びます:

理論

実践

レファレンス


理論

クラスタリングと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

注:化合物は次の場合に同じクラスターに含められます。(距離行列を使用する場合)クラスターのセントロイドからの最長距離が特定のカットオフ値以下の時、あるいは、(類似度行列を使用する場合)最小の類似度が特定のカットオフ値以上のときです。

  • 潜在的クラスター中心の発見

    • クラスター中心は、与えられたクラスターの中で最も近傍化合物の数が多い分子です。
    • 近傍化合物の割り当て(Annotate neighbors):各分子について、与えられた閾値以下のタニモト距離に含まれる全ての分子の数を数えます。
    • 化合物を近傍化合物の数によって降順に並べ替え、クラスター中心となりうるもの(即ち、近傍の数が最大となる化合物)をファイルの一番上に置くようにします。
  • 排除球(exclusion shpheres)に基づくクラスタリング

    • 並べ替えたリスト内の最初の分子(中心、centroid)から始めます。
      • クラスタリングに使われたカットオフ値以上のタニモトインデックスをもつ全ての化合物をクラスターのメンバーとします(類似度の場合)
        • 考慮しているクラスターのメンバーとなった各化合物にフラグをたて、以降の比較の対象から取り除きます。したがって、フラグがたてたれた化合物は他のクラスターの中心やメンバーとなることはできません。この手順は新しく作成したクラスターの周りに排除球を置くようなものです。
        • リストの最初の化合物について全ての近傍化合物がみつかったら、リストトップの最初に利用可能な(即ちフラグが立っていない)化合物を新しいクラスター中心とします。
      • リストのフラグの立っていない全ての化合物について、上から順に同じ手順を繰り返していきます。
    • クラスタリング手順の最後までフラグの立てられなかった化合物はシングルトン(singleton)となります。
      • シングルトンとして割り当てられた化合物の中には、指定したタニモト類似度インデックスで近傍となる化合物をもつものありえますが、これらの近傍化合物は、クラスター中心に基づく枠組みによってすでに排除球に含められていることに注意してください。
from IPython.display import IFrame
IFrame('./images/butina_full.pdf', width=600, height=300)

f:id:magattaca:20200425151951p:plain 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()

f:id:magattaca:20200425152732p:plain

合理的なカットオフ値を選ぶにはどうすればよいか?

クラスタリングの結果はユーザーの選んだ閾値に依存するので、カットオフ値の選択についてより詳細に見てみたいと思います。

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()

f:id:magattaca:20200425152815p:plain

f:id:magattaca:20200425152827p:plain

f:id:magattaca:20200425152837p:plain

f:id:magattaca:20200425152846p:plain

f:id:magattaca:20200425152857p:plain

f:id:magattaca:20200425152910p:plain

f:id:magattaca:20200425152922p:plain

f:id:magattaca:20200425152933p:plain

f:id:magattaca:20200425152943p:plain

f:id:magattaca:20200425152952p:plain

結果を見ると、閾値(距離のカットオフ)が高いほど、より多くの化合物が類似しているとみなされ、したがってより少ないクラスター数にクラスタリングされます。閾値が低くなると、より多数の小さなクラスターとシングルトン(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

f:id:magattaca:20200425153106p:plain

クラスターの可視化

最大のクラスターからの化合物例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:

f:id:magattaca:20200425153234p:plain

# トークトリアル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:

f:id:magattaca:20200425153306p:plain それぞれのクラスターの最初の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:

f:id:magattaca:20200425153341p:plain

最初の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')
# 平均=赤色, 中央値=青色

f:id:magattaca:20200425153421p:plain

化合物選択

以下では、多様性のあるサブセットとして最大 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()

f:id:magattaca:20200425153608p:plain

クイズ

訳以上

追記

今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。

kNNによるクラスタリング:
* 化学の新しいカタチ:RDKitとscikit-learnで機械学習:変異原性をk-最近傍法で予測

Butinaクラスタリングは以下でも取り上げられています:
* AI創薬のためのケモインフォマティクス 入門:6章化合物の類似性を評価してみる

多様なライブラリの構築とその方法については以下で説明してくださっています:
* 化学の新しいカタチ:ケモインフォマティクスで多様な化合物ライブラリーを構築する

また、今回のトークトリアルでは取り上げられていませんでしたが、クラスタリングと関連する可視化として次元圧縮によるプロットがあげられると思います。ケミカルスペースの可視化については以下で解説してくださっています。

勝手にいろいろと引用させていただいておりますので問題があれば削除します。

TeachOpenCADD トピック4 〜リガンドベーススクリーニング〜

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

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

以下、訳。

トークトリアル 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類似性
  • 類似性検索によるバーチャルスクリーニング
    • データセットの全化合物に対する化合物クエリの比較
    • 類似度の分布
    • 最も類似した分子の描画
    • エンリッチメントプロットの作成

レファレンス


理論

化合物類似性

化合物類似性は化学情報学(ケモインフォマティクス、chemical informatics)の中でよく知られており、頻繁に用いられる考え方です。化合物とその特性の比較はいろいろな用途で使用でき、望みの特性と生理活性をもつ新しい化合物を見つけるのに役立つかもしれません。

構造的に類似の化合物は類似の特性、そして類似の生理活性を示すという考え方は、類似性質原則(similar property principle、SPP)や構造活性相関(structure activity relationship、SAR)に表れています。この文脈において、バーチャルスクリーニングは、結合親和性のわかっている化合物セットがあれば、そのような化合物をさらに探すことができる、というアイデアに基づいています。

化合物記述子

類似度は適用範囲に応じて様々な方法で評価することができます(参考 J. Med. Chem. (2014), 57, 3186-3204):

  • 1D 化合物記述子: 溶解度、logP、分子量、融点 etc.
    • グローバル記述子(Global descriptor):分子全体を一つの値だけで表現する
    • 通常、機械学習(machine learning、ML)を適用するには分子を特定するのに十分な特性とはならない
    • 機械学習のための化合物エンコーディングを改良するために2Dフィンガープリントに付け加えることができる
  • 2D 化合物記述子: 分子グラフ(Molecular graph)、経路(path)、フラグメント、原子環境(atom environment)
    • 分子の個々の部位の詳細な表現
    • 一つの分子に対して多数のフィンガープリントと呼ばれる特徴/ビット
    • 類似性検索と機械学習で非常によく使われる
  • 3D 化合物記述子: 形状(Shape), 立体化学
    • 化学者は通常2次元表現で訓練されている 
    • 化合物の自由度(flexibility、化合物の「正しい」配座はどれか?)のため、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).

f:id:magattaca:20200419214525p:plain:w350
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)を参照してください。

f:id:magattaca:20200419214559p:plain:w350
Figure 3: Morganサーキュラーフィンガープリントの図例(Andrea Morgerによる図)

化合物類似性評価

記述子/フィンガープリントの計算ができれば、それらを比較することで、二つの分子の間の類似度を評価することができます。化合物類似度は様々な類似度係数で定量化することができますが、よく使われる2つの指標はタニモト係数とDice係数です(Tanimoto and Dice index) (J. Med. Chem. (2014), 57, 3186-3204)。

タニモト係数(Tanimoto coefficient)


T _{c}(A,B) = \frac{c}{a+b-c}

a: 化合物Aに存在する特徴の数 
b: 化合物Bに存在する特徴の数 
c: 化合物AとBで共有されている特徴の数

Dice係数(Dice coefficient)


D_{c}(A,B) = \frac{c}{\frac{1}{2}(a+b)}

a: 化合物Aに存在する特徴の数 
b: 化合物Bに存在する特徴の数 
c: 化合物AとBで共有されている特徴の数

類似度評価は通常、それぞれのフィンガープリントの正の(1の)ビットの数と、両者が共通してもつ正のビットの数を考慮します。Dice類似度は通常タニモト類似度よりも大きな値を返し、それはそれぞれの分母の違いに由来します。:


\frac{c}{a+b-c} \leq \frac{c}{\frac{1}{2}(a+b)}

バーチャルスクリーニング(Virtual screening)

医薬品探索の初期段階における課題は、低分子(化合物)のセットを、有りうる巨大なケミカルスペースから、研究対象のターゲット分子に結合するポテンシャルのあるものに範囲を狭めることです。このケミカルスペースは非常に大きく、低分子化合物群は部分構造(chemical moiety)の1020 の組み合わせにまでいたります (ACS Chem. Neurosci. (2012), 19, 649-57) 。

目的のターゲット分子に対するこれら低分子の活性を評価するハイスループットスクリーニング(HTS)実験は費用と時間が非常にかかるので、計算機に支援された(computer-aided)手法により、試験にかける低分子のリストをより絞り込む(focused list)ことが期待されています。このプロセスはバーチャル(ハイスループット)スクリーニングと呼ばれていて、研究対象のターゲット分子に結合する可能性の最も高い低分子を見つけるために、巨大な低分子ライブラリーをルールとパターンのどちらか、あるいは両方によってフィルタリングします。

類似度検索を用いたバーチャルスクリーニング

バーチャルスクリーニングの簡単な方法として、既知の活性化合物(群)と新規な化合物セットを比較して、最も類似しているものを探すことが行われます。類似性質原則(similar property principle、SPP)に基づき、(例えば既知の阻害剤に)最も類似した化合物は類似の効果を有すると推測されます。類似性検索に必要となるものは次の通りです(上でも詳細に議論しています)。

  • 化学/分子の特徴をエンコードした表現
  • 特徴のポテンシャルの重み付け(オプション)
  • 類似度評価

類似性検索はある特定のデータベースの全ての化合物と一つの化合物との間の類似度を計算することで実行することができます。データベースの化合物の類似度係数によるランク付けにより、最も類似度の高い分子が得られます。

エンリッチメントプロット(Enrichment plots)

エンリッチメントプロットはバーチャルスクリーニングの結果の妥当性を評価するために使われ、ランク付けされたリストの上位x%の中に含まれる活性化合物の比率を表します。すなわち、

  • データセット全体のうち、トップにランクした化合物の比率(x-axis)  vs.
  • データセット全体のうち活性化合物(y-axis)の比率

[f:id:magattaca:20200419214803p:plain: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)

f:id:magattaca:20200419215200p:plain

化合物記述子の計算

化合物の比較を行うために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))

f:id:magattaca:20200419215245p:plain

見てわかるように類似の分子量を持つ化合物は類似の構造をもつことがあります(例  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つの類似度評価、すなわちTanimotoDiceを、2つのタイプのフィンガープリント、すなわちMACCSMorganフィンガープリントに適用します。

例: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")

f:id:magattaca:20200419215752p:plain

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")

f:id:magattaca:20200419215839p:plain

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()

f:id:magattaca:20200419215900p:plain

異なるフィンガープリント(ここでは、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

f:id:magattaca:20200419215930p:plain

# クエリ化合物の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()

f:id:magattaca:20200419220014p:plain

ここでも類似度を比較します。今回は直接、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()

f:id:magattaca:20200419220032p:plain

類似度分布は類似度を解釈するのに重要です(例 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)

f:id:magattaca:20200419220106p:plain

データセットにおいてゲフィチニブと比較してトップにランクした化合物群は、最初は我々のデータセットに含まれるゲフィチニブのエントリー(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()

f:id:magattaca:20200419220159p:plain

エンリッチメントプロットによると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)となっています。この値は先のRandomOptimalと比較して妥当な値に思います。

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の比較

私は既存のメソッドの使い方ですらよくわからん、誰か日本語で解説して!という感じなので実際にいろいろ作ってしまう方々は本当にすごいなあと思います。アホですみません。

勝手にいろいろと引用させていただいておりますので問題があれば削除します。