magattacaのブログ

日付以外誤報

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

TeachOpenCADD

こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル４をもととしており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）を適用するには分子を特定するのに十分な特性とはならない
    • 機械学習のための化合物エンコーディングを改良するために２Dフィンガープリントに付け加えることができる
  • 2D 化合物記述子: 分子グラフ（Molecular graph）、経路（path）、フラグメント、原子環境（atom environment）
    • 分子の個々の部位の詳細な表現
    • 一つの分子に対して多数のフィンガープリントと呼ばれる特徴/ビット
    • 類似性検索と機械学習で非常によく使われる
  • 3D 化合物記述子: 形状（Shape）, 立体化学
    • 化学者は通常２次元表現で訓練されている　
    • 化合物の自由度（flexibility、化合物の「正しい」配座はどれか？）のため、２次元表現と比べて頑健性が低い
  • 生物学的類似性
    • 生物学的フィンガープリント（例、個々のビットが異なるターゲット分子に対して評価された生理活性を表す）
    • 化合物構造からは独立
    • 実験データ（あるいは予測値）が必要

すでに トークトリアル　2　で、分子量やlogPといった1D　物理化学パラメーターを計算する方法を学びました。RDKitに実装されているそのような記述子は RDKit documentation: Descriptors　で見つけることができます。

以降では、２D（あるいは３D）化合物記述子の定義に焦点を当てます。多くの場合、分子ごとに固有の（ユニークな）ものとなるので、これらの記述子はフィンガープリント（指紋）ともよばれます。

化合物フィンガープリント（Molecular fingerprints）

部分構造に基づくフィンガープリント（Substructure-based fingerprints）

化合物フィンガープリントは化学的な特徴と分子の特徴をビット文字列（bitstring）やビットベクトル（bitvector）、配列（array）の形でエンコードします。各ビットは、事前に定義された分子の特徴あるいは環境に相当し、「１」はその特徴が存在していることを、「０」は存在していないことを示します。実装の仕方によっては、数え上げベース（count-baseｄ）となっていて、ある特定の特徴がいくつ存在しているかを数えるようになっていることに注意してください。

フィンガープリントのデザインには複数の方法があります。ここではよく使われる２Dフィンガープリントのとして、MACCSキーとMorganフィンガープリントの２種類を導入します。 RDKit documentation: Fingerprintsに記載されているように、RDKitではこの２つ以外にも多数のフィンガープリントを提供しています。

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による図)

化合物類似性評価

記述子/フィンガープリントの計算ができれば、それらを比較することで、二つの分子の間の類似度を評価することができます。化合物類似度は様々な類似度係数で定量化することができますが、よく使われる２つの指標はタニモト係数と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で共有されている特徴の数

類似度評価は通常、それぞれのフィンガープリントの正の（１の）ビットの数と、両者が共通してもつ正のビットの数を考慮します。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を使って化合物のエンコード（化合物フィンガープリント）をしたのち、上の理論編で議論したように、類似度（化合物類似性評価）を計算するため、それらの比較を実行します。

２番目のパートではこれらのエンコーディングと比較手法を使って類似度検索（バーチャルスクリーニング）を実施します。既知のEGFR阻害剤ゲフィチニブ（Gefitinib)をクエリとして使用し、EGFRに対して試験済みの化合物データセットの中から類似した化合物を検索します。このデータセットは トークトリアル１でChEMBLから抽出し、トークトリアル2でリピンスキーのルールオブファイブによりフィルタリングをおこなったものです。

化合物の読み込みと描画

まず、８個の化合物例を定義し描画します。後ほど、これらの分子をエンコードし比較します。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と２Dの化合物記述子を抽出、生成します。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)。

次に、より詳細な分子の特徴を説明するために、２D化合物記述子を見て見ましょう。

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サーキュラーフィンガープリントも計算します。２つの異なる関数により、Morganフィンガープリントは整数（int）あるいはビット（bit）ベクトルとして計算することができます。

# Morganフィンガープリントを生成（int vector）、デフォルトでは半径２でベクトルの長さは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として）生成、デフォルトでは半径２でフィンガープリントの長さは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)

化合物類似度の計算

次では、２つの類似度評価、すなわちTanimotoとDiceを、２つのタイプのフィンガープリント、すなわちMACCSとMorganフィンガープリントに適用します。

例：２つのMACCSフィンガープリントをタニモト類似度で比較

# ２つの化合物のタニモト係数を計算
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)で、次の分子はランク１から始まる
    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を最も似ていない化合物として認識しました。一方で、２番目にランク付されたのはHemi-cycline Dでした。この化合物はサイクリン系化合物の部分構造で、Morganフィンガープリントのアルゴリズムが原子の環境に基づくものであることがその理由であるかもしれません（一方で、MACCSフィンガープリントは特定の特徴の出現頻度を求めるものとなっています）。

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

フィンガープリントと類似度の計算方法を学んだので、この知識を化合物セット全体からのクエリ化合物の類似度検索に応用することができます。

既知のEGFR阻害剤ゲフィチニブ（Gefitinib)をクエリとして使用し、EGFRに対して試験済みの化合物データセットの中から類似した化合物を検索します。このデータセットは トークトリアル１でChEMBLから抽出し、トークトリアル2でリピンスキーのルールオブファイブによりフィルタリングをおこなったものです。

クエリ化合物をデータセットの全化合物と比較

トークトリアル2で取得したChEMBLデータベースから取り出したEGFRに対して評価済みのフィルタリングされた化合物を含むcsvファイルから化合物を読み込みます。１つのクエリ化合物（ここではゲフィチニブ）を使って、類似の化合物をデータセットの中から探し出します。

# 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類似度の値よりも小さくなります。また、２つの異なるフィンガープリント（例　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

ここでも類似度を比較します。今回は直接、２つのフィンガープリントに関するタニモト類似度と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は２つの類似度評価基準のデータを繋げたものなので、それぞれ別々にしてみます。

# tanimoto_MACCS
enrich_df_taMA = enrich_df[enrich_df['similarity_measure'] == 'tanimoto_MACCS']

# tanimoto_morgan
enrich_df_tamo = enrich_df[enrich_df['similarity_measure'] == 'tanimoto_morgan']

print("Size of enrich_df: ", len(enrich_df))
print("Size of tanimoto_MACCS DataFrame: ", len(enrich_df_taMA))
print("Size of tanimoto_morgan DataFrame: ", len(enrich_df_tamo))

# Size of enrich_df:  9046
# Size of tanimoto_MACCS DataFrame:  4523
# Size of tanimoto_morgan DataFrame:  4523

DataFrameの5% ranked Datasetに相当する箇所を見てみます。

# 5% に相当する数
index_5perc = round(len(enrich_df_taMA)*0.05)

# DataFrameのindexは0から始まるので-1した行を表示
enrich_df_taMA[index_5perc-1:index_5perc]
% ranked dataset % true actives identified similarity_measure
225 0.049967 0.07319 tanimoto_MACCS

見やすさのためsliceでDataFrameとして取り出しています。
ランク上位5%(0.049)に相当する数のなかに、実際の活性評価でactiveだったものは7.3%(% true actives identified, 0.07319)となっています。この値は先のRandom、Optimalと比較して妥当な値に思います。

DataFrameのデータ自体には問題なさそうなので、値の取り出し方（関数print_data_ef）に問題があったのでしょう。 関数の中身を順番に実行してみます。

# 5%に設定
perc_ranked_dataset = 5

# 5%内のDataFrameを取り出し、その一番最後の行(tail)を取り出す
enrich_df[enrich_df['% ranked dataset'] <= perc_ranked_dataset].tail(1)
% ranked dataset % true actives identified similarity_measure
4522 1.0 1.0 tanimoto_morgan

取り出されたのは「index:4522」で4523番目の行です。この% true actives identifed列がEFとして取り出されていた値(1.0)です。
閾値5以下で取り出されたのはsimilarity_measure:tanimoto_morganの全化合物でした。単純にDataFrameのデータは%に換算していないのに、取り出す際に%換算の値を使ってしまったのが原因のようです。

それではそれぞれの類似度について正しい値を確認してみます。

# 関数の再定義
def print_data_ef2(perc_ranked_dataset, enrich_df):
    perc_ranked_dataset_100 = perc_ranked_dataset / 100
    data_ef = enrich_df[enrich_df['% ranked dataset'] <= perc_ranked_dataset_100].tail(1)
    data_ef = round(float(data_ef['% true actives identified'] * 100), 1)
    print("Experimental EF for ", perc_ranked_dataset, "% of ranked dataset: ", data_ef, "%", sep="")

# MACCS keyの場合
# パーセンテージを選択
perc_ranked_list = 5

# EFデータを取得
print_data_ef2(perc_ranked_list, enrich_df_taMA)
print_random_ef(perc_ranked_list)
print_optimal_ef(perc_ranked_list, similarity_df, threshold)

# Experimental EF for 5% of ranked dataset: 7.3%
# Random EF for 5% of ranked dataset:       5.0%
# Optimal EF for 5% of ranked dataset:      8.9%

# Morganフィンガープリントの場合
# パーセンテージを選択
perc_ranked_list = 5

# EFデータを取得
print_data_ef2(perc_ranked_list, enrich_df_tamo)
print_random_ef(perc_ranked_list)
print_optimal_ef(perc_ranked_list, similarity_df, threshold)

#  Experimental EF for 5% of ranked dataset: 7.9%
#  Random EF for 5% of ranked dataset:       5.0%
#  Optimal EF for 5% of ranked dataset:      8.9%

いずれも「optimal > Experimental > Random」となっており、Morganの方がMACCSよりも若干良い値となっています。無事エンリッチメントプロットと比較してもおかしくない値が得られました。

訳注ここまで

ディスカッション

ここではタニモト類似度を使ってバーチャルスクリーニングを実施しました。もちろん、Dice類似度や他の類似度評価を使っても行うことができます。

化合物フィンガープリントを使用した類似度検索の欠点は、化合物類似度に基づくものなので新規な構造を生み出さないことです。化合物類似度を扱う上でのその他の課題としては、いわゆるアクティビティクリフ（activity clliff）があります。分子の官能基におけるちょっとした変化が生理活性の大きな変化を起こすことがあります。

クイズ

  • アクティビティクリフを回避するにはどこから始めれば良いでしょうか？
  • MACCSフィンガープリントとMorganフィンガープリントを互いに比較した場合の利点と欠点は何でしょう？
  • 使用したフィンガープリントによっておこる、類似度データフレームにおける順序の違いをどう説明できるでしょうか？

訳以上

追記

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

フィンガープリントと分子の類似性について:
* 化学の新しいカタチ：RDKitでフィンガープリントを使った分子類似性の判定

また、同書の7章: グラフ構造を利用した類似性の評価にはActivity Cliffの説明と事例も取り上げられています。全体を通して非常に勉強になるので未読の方は是非一読をお勧めします。

エンリッチメントファクターとその他の指標については以下の記事が参考になりました。
tonetsの日記: バーチャルスクリーニングの指標のテスト その2

難しすぎて私では理解できませんでしたがECFPを実装されている方もいらっしゃいました。
創薬・材料探索のための機械学習: ECFPとNeural-fingerprintの比較

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

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

TeachOpenCADD トピック３ 〜部分構造による化合物フィルタリング〜

TeachOpenCADD

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

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

github.com

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

以下、訳。

トークトリアル 3

化合物フィルタリング：好ましくない部分構造

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

Maximilian Driller and Sandra Krüger

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

いくつか私たちのスクリーニングライブラリーに含めたくない部分構造があります。このトークトリアルでは、そのような好ましくない部分構造の様々なタイプを学び、そしてRDKitを使ってそれらの部分構造を見つけ、ハイライトする方法を学びます。

学習の目標

理論

  • 好ましくない部分構造とは何か?
  • 汎分析干渉化合物 Pan Assay Interference Compounds (PAINS)

実践

  • ChEMBLデータベースから化合物のセットを読み込む（トークトリアル 2　で準備したもの）
  • RDKitでの実装を利用し好ましくない部分構造をフィルターにかけ取り除く
  • 独自の好ましくない部分構造のリストを作成し、フィルタリングを実施
  • 部分構造の検索とハイライト

レファレンス

理論

好ましくない部分構造

部分構造には好ましくないものがあります。例えば毒性あるいは反応性があるといった理由や、薬物動態学的特性が好ましくないという理由、あるいは特定のアッセイに干渉する可能性が高いという理由です。 今日の医薬品探索ではよくハイスループットスクリーニング (HTS wikipedia)　を実施します。好ましくない部分構造をフィルタリングすることで、より望ましいスクリーニングライブラリを構築することができます。これにより、スクリーニングの前にライブラリの数を減らすことができ、時間と資源の節約につながります。

Brenkらは顧みられない病気の治療のための化合物の探索に使用するスクリーニングライブラリーをフィルタリングするため、好ましくない部分構造のリストを作成しました(Chem. Med. Chem. (2008), 3,435-444)。好ましくない特徴の例として、ニトロ基（変異原性の問題）、スルホン酸基とリン酸基（薬物動態学的特性が好ましくない可能性が高い）、2-ハロピリジンとチオール基（反応性の問題）、といったものが挙げられます。

好ましくない部分構造のリストは上記の文献に報告されており、このトークトリアルの実践編でも使用します。

汎分析干渉化合物 Pan Assay Interference Compounds (PAINS)

概要

PAINS (PAINS wikipedia) は実際には偽陽性にも関わらず、HTSでしばしばヒットとして見出される化合物です。PAINSは特定の一つのターゲット分子と反応するというよりもむしろ、非特異的に多数のターゲット分子と反応する傾向があります。通常、非特異的な結合もしくはアッセイの構成要素との相互作用により、様々なアッセイで様々なタンパク質に対して見かけ上結合します。

f:id:magattaca:20200418212342j:plain

Figure 1:汎分析干渉化合物（Pan Assay Interference Compounds、PAINS）の観点における特異的結合と非特異的結合。図はwikipediaより引用

Baellらによる文献で用いられたフィルター

Baellらはアッセイのシグナルに干渉する部分構造に焦点を当てました (J. Med. Chem. (2010), 53(7),2719-2740)。そのようなPAINSを見つけるのに役立つ部分構造について記述し、部分構造フィルタリングに利用できるリストを提供しました。

実践

データの読み込みと可視化

まず、必要なライブラリをインポートし、トークトリアル T2 でフィルタリングしたデータセットを読み込み、はじめの分子群を描画します。

import pandas
from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import rdFMCS
from rdkit.Chem import AllChem
from rdkit.Chem import Descriptors
from rdkit.Chem import Draw
from rdkit import DataStructs
from rdkit.Chem import PandasTools
import matplotlib.pyplot as plt

filteredData = pandas.read_csv("../data/T2/EGFR_compounds_lipinski.csv", delimiter=";", index_col=0)
filteredData.drop(['HBD','HBA','MW','LogP'], inplace=True, axis=1) # 不必要な情報の除去
print ('Dataframe shape: ', filteredData.shape) # データフレームの次元をプリント
filteredData.head(5)

# Dataframe shape:  (4523, 6)
molecule_chembl_id units IC50 smiles pIC50 rule_of_five_conform
0 CHEMBL63786 nM 0.003 Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 11.522879 yes
1 CHEMBL53711 nM 0.006 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849 yes
2 CHEMBL35820 nM 0.006 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849 yes
3 CHEMBL53753 nM 0.008 CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.096910 yes
4 CHEMBL66031 nM 0.008 Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 11.096910 yes 
PandasTools.AddMoleculeColumnToFrame(filteredData, smilesCol='smiles') # 分子を列に追加
# 最初の２０個の分子を描画
Draw.MolsToGridImage(list(filteredData.ROMol[0:20]), 
                    legends=list(filteredData.molecule_chembl_id[0:20]), 
                    molsPerRow=4)

f:id:magattaca:20200418212826p:plain

RDKitを使用してPAINSをフィルタリング

PAINSフィルターはすでにRDKitに実装されています(RDKit Documentation)。使ってみましょう！

from rdkit.Chem.FilterCatalog import *
params = FilterCatalogParams()
# 全　PAINS (A, B and C)　からカタログを構築
params.AddCatalog(FilterCatalogParams.FilterCatalogs.PAINS)
catalog = FilterCatalog(params)

# フィルタリングしたデータを格納するための空のデータフレームを作成
rdkit_highLightFramePAINS = pandas.DataFrame(columns=('CompID', 'CompMol', 'unwantedID'))
rdkit_noPAINS = pandas.DataFrame(columns=('ChEMBL_ID', 'smiles','pIC50'))
rdkit_withPAINS = pandas.DataFrame(columns=('ChEMBL_ID', 'smiles', 'pIC50','unwantedID'))

# フィルタリングしたデータフレームのインデックスと行についてforループを回す
for i,row in filteredData.iterrows():
    curMol = Chem.MolFromSmiles(row.smiles) # 現在の分子
    match = False # Falseにmatchを設定
    rdkit_PAINSList = []
    # 最初のmatchを取得
    entry = catalog.GetFirstMatch(curMol)
    
    if entry!=None:
        # 現在の好ましくない部分構造の名前をリストに追加
        rdkit_PAINSList.append(entry.GetDescription().capitalize())
        # データフレームに関連するマッチング情報を追加
        rdkit_highLightFramePAINS.loc[len(rdkit_highLightFramePAINS)] = [row.molecule_chembl_id, curMol,
        entry.GetDescription().capitalize()]
        match = True
    if not match:
        # PAINSを含まない化合物のデータフレームに追加
        rdkit_noPAINS.loc[len(rdkit_noPAINS)] = [row.molecule_chembl_id, row.smiles, row.pIC50]
    else: 
        # PAINSを含む化合物のデータフレームに追加
        # 好ましくない部分構造を含むデータフレームに関連する情報を加える
        rdkit_withPAINS.loc[len(rdkit_withPAINS)] = [row.molecule_chembl_id, row.smiles, row.pIC50, entry.GetDescription().capitalize()]

print("Number of compounds with PAINS: %i"%(len(rdkit_withPAINS)))
print("Number of compounds without PAINS: %i (=remaining compounds)"%(len(rdkit_noPAINS)))

rdkit_highLightFramePAINS.head(10)

# Number of compounds with PAINS: 394
# Number of compounds without PAINS: 4129 (=remaining compounds)
CompID CompMol unwantedID
0 CHEMBL93032 Mol Catechol_a(92)
1 CHEMBL2029429 Mol Anil_di_alk_a(478)
2 CHEMBL2029428 Mol Anil_di_alk_a(478)
3 CHEMBL4071474 Mol Anil_di_alk_a(478)
4 CHEMBL4126630 Mol Anil_di_alk_a(478)
5 CHEMBL3663926 Mol Anil_di_alk_a(478)
6 CHEMBL3663925 Mol Anil_di_alk_a(478)
7 CHEMBL3663936 Mol Anil_di_alk_a(478)
8 CHEMBL3663923 Mol Anil_di_alk_a(478)
9 CHEMBL3663934 Mol Anil_di_alk_a(478)

訳注(04/2020)

訳者はPythonに詳しくないのでrdkit_highLightFramePAINS.loc[len(rdkit_highLightFramePAINS)] の部分に関する備忘録です。
locはPandas DataFrameで行をindexで指定する方法です。 len(~)でDataFrameの行数を取得することですでにデータが入っている行数がわかります。
DataFrameのindexは0から始まっているので、DataFrameの行数をindexとして指定することで、新しいデータを新しい行（既存の行数+1の行）に追加することになるのだと思います。

また、PAINSの利用例はRDKit Documentationのこちらに同様のコードが見られます。

訳注ここまで

好ましくない/毒性の懸念される部分構造でフィルタリング(Brenkのリスト)

RDKitには、PAINSのような好ましくない部分構造のリストがいくつかすでに実装されています。ですが、実装されていないRDKit外部のリストを使って部分構造の一致検索を行うことも可能です。ここではBrenkらによる文献(Chem. Med. Chem. (2008), 3,435-444)のsupporting informationで提供されているリストを使用します。

注）　データをダウンロードし、dataフォルダにcsvファイルとして保存済みです（形式は name-space-SMARTS です）。まず、データを読み込みます。

unwantedSubs = []
unwantedNames = []

for line in open('../data/T3/unwantedSubstructures.csv', 'r'):
    if not line.startswith("#"): # ヘッダーを無視
        splitted = line.strip().split(" ") # 各行を分割
        m = Chem.MolFromSmarts(splitted[1]) # SMARTSから分子を生成
        name = splitted[0].capitalize() # nameに名前を保存
        unwantedNames.append(name) # 好ましくない部分構造の名前をリストに追加
        unwantedSubs.append(m) # 好ましくない部分構造をリストに追加
print("Number of unwanted substructures in list =", len(unwantedSubs)) # 好ましくない部分構造の数を表示

# Number of unwanted substructures in list = 104

訳注(04/2020)

Pythonやプログラミング初心者仲間のための蛇足・・・
「unwantedSubstructures.csv」ファイルは１行目のみ「#name smart」となっています。if not line.starwith(#)の部分で「#」から始まる１行目(header)を飛ばしています。
strip()は文字列の先頭・末尾の余分な文字を削除するメソッドで、空白文字としてスペース、タブ、改行も取り除けるのでcsvの改行文字を除いているのだと思います。
split(" ")はスペースで分割しているので、改行文字を取り除いた各行のname-space-SMARTS形式をspaceで分割しnameとSMARTSにしているのだと思います。
こう考えるとsplitted[0]がname、splitted[1]がSMARTSになっていることがわかりやすい気がします。

訳注ここまで

２、３個部分構造を見てみましょう　(すべてのSMARTSが表示できるわけではありません。一部を選んでいます。)。

Chem.Draw.MolsToGridImage(list(unwantedSubs[2:5]), subImgSize=(200, 300), legends=unwantedNames[2:5])

f:id:magattaca:20200418215706p:plain

これらの好ましくない部分構造とマッチするものがあるか、フィルタリングしたデータフレームの中を検索してみます。

# フィルタリングしたデータを格納するためのデータフレームを作成
highLightFrameUNW = pandas.DataFrame(columns=('CompID', 'CompMol', 'unwantedID', 'unwSubstr'))
noUnwanted = pandas.DataFrame(columns=('ChEMBL_ID', 'smiles','pIC50'))
withUnwanted = pandas.DataFrame(columns=('ChEMBL_ID', 'smiles', 'pIC50','unwantedID'))
molsToDraw = []

# データセットの各化合物に対して
for i,row in filteredData.iterrows(): # フィルタリングしたデータフレームのインデックスと行についてfor ループを回す
    curMol = Chem.MolFromSmiles(row.smiles) # 現在の分子
    match = False # Falseにmatchを設定
    unwantedList = []
    molsToDraw.append(curMol)
    
    # 全ての好ましくない部分構造を検索
    for idx, unwSub in enumerate(unwantedSubs):
        # 部分構造があるかチェック
        if curMol.HasSubstructMatch(unwSub): # 現在の分子が好ましくない部分構造を有する場合
            match = True # matchをTrueに設定
            unwantedList.append(unwantedNames[idx]) # 好ましくない部分構造の名前をリストに追加
            # 関連する情報をデータフレームに追加
            highLightFrameUNW.loc[len(highLightFrameUNW)] = [row.molecule_chembl_id, curMol, unwantedNames[idx], unwSub]
    if not match: # 一致する部分構造が見つからない場合
        noUnwanted.loc[len(noUnwanted)] = [row.molecule_chembl_id, row.smiles, row.pIC50]
        # 欲しい部分構造のデータフレームに関連情報を追加
    else: # 一致が見つかった場合
        withUnwanted.loc[len(withUnwanted)] = [row.molecule_chembl_id, row.smiles, row.pIC50, unwantedList] #　好ましくない部分構造データフレームに関連する情報を追加

print("Number of compounds with unwanted substructures: %i"%(len(withUnwanted)))
print("Number of compounds without unwanted substructures: %i (=remaining compounds)"%(len(noUnwanted)))

# Number of compounds with unwanted substructures: 2356
# Number of compounds without unwanted substructures: 2167 (=remaining compounds)

highLightFrameUNW.head(4)
CompID CompMol unwantedID unwSubstr
0 CHEMBL63786 Mol Polycyclic-aromatic-hydrocarbon Mol
1 CHEMBL1243316 Mol Aliphatic-long-chain Mol
2 CHEMBL1243316 Mol Michael-acceptor Mol
3 CHEMBL1243316 Mol Triple-bond Mol

部分構造を分子の中で直接ハイライトすることもできます。

first_highLightFrameUNW = highLightFrameUNW.head(8) # リストの最初の8エントリーのサブセット

# 分子を描画し、好ましくない部分構造をハイライトする 
Draw.MolsToGridImage(list(first_highLightFrameUNW["CompMol"]), subImgSize=(400,300),
    molsPerRow=2, highlightAtomLists=
    [m.GetSubstructMatch(first_highLightFrameUNW["unwSubstr"][i]) for i,m in enumerate(first_highLightFrameUNW["CompMol"])], 
    legends=list(first_highLightFrameUNW["CompID"]+": "+first_highLightFrameUNW["unwantedID"]))

f:id:magattaca:20200418215814p:plain

訳注(04/2020)

Pythonやプログラミング初心者仲間のための蛇足・・・
[m.GetSubstructMatch(first_highLightFrameUNW["unwSubstr"][i]) for i,m in enumerate(first_highLightFrameUNW["CompMol"])]の部分ですが、enumerateはインデックス番号と要素を同時に取得するメソッドで、ここではそれぞれiとmを割り当てています。
要素m分子がもつ好ましくない部分構造を見つけるために、RDKitのGetSubstructMatchを使用しており、検索対象の好ましくない部分構造をunwSubstr列のインデックスiを参照するという形で実行されています。

訳注ここまで

例をSVGファイルとして保存します。

# イメージをファイルに保存
img = Draw.MolsToGridImage(list(first_highLightFrameUNW["CompMol"]), subImgSize=(400,300),
    molsPerRow=3, highlightAtomLists=
    [m.GetSubstructMatch(first_highLightFrameUNW["unwSubstr"][i]) for i,m in enumerate(first_highLightFrameUNW["CompMol"])], 
    legends=list(first_highLightFrameUNW["unwantedID"]), useSVG=True)

# SVGデータの取得
molsvg = img.data

# 不透明な背景を透明に置き換える
molsvg = molsvg.replace("opacity:1.0", "opacity:0.0");
molsvg = molsvg.replace("12px", "24px");

# 変更されたSVGデータをファイルに保存
f = open("../data/T3/substructures.svg", "w")
f.write(molsvg)
f.close()

好ましくない部分構造を持つ化合物と持たない化合物のリストを保存。

# 好ましくない部分構造を有する化合物をcsvファイルに書き出し
withUnwanted.to_csv("../data/T3/EGFR_compounds_lipinski_noPAINS.csv", sep=',') 

# 好ましくない部分構造を持たない化合物をcsvファイルに書き出し
noUnwanted.to_csv("../data/T3/EGFR_compounds_lipinski_noPAINS_noBrenk.csv", sep=',') 

# 好ましくない部分構造を有する化合物のcsvファイルの最初のいくつかを表示
noUnwanted.head()
ChEMBL_ID smiles pIC50
0 CHEMBL53711 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849
1 CHEMBL35820 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849
2 CHEMBL53753 CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.096910
3 CHEMBL66031 Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 11.096910
4 CHEMBL176582 Cn1cnc2cc3ncnc(Nc4cccc(Br)c4)c3cc21 11.000000

見つかった好ましくない部分構造をさらに解析します。

# 最も頻度の高い化合物の数を数えます。
unwCounts = {}
for ele in highLightFrameUNW.values:
    unw = ele[2] # highLightFrameUNWデータフレームに含まれる好ましくない部分構造のID
    if unwCounts.get(unw, "empty") == "empty": # 好ましくない構造のIDがまだ辞書に含まれていない場合
        unwCounts[unw] = [1, ele[3]] # 1を付与し、辞書の好ましくない構造に追加
    else: # もしkey （好ましくない構造のID）がすでに存在しているなら、出現回数の値を１増やす 
        unwCounts[unw][0] += 1

frequentUNW = []
frequentUNWNames = []

# unwCountsに含まれる構造: IDをkey、出現頻度と分子をvalueとしてもつ辞書
# 例　('acyclic-C=C-O', [7, <rdkit.Chem.rdchem.Mol object at 0x7fa58fc06710>])

# 部分構造の頻度で辞書をソートする
for key, value in sorted(unwCounts.items(), key=lambda kv: kv[1][0], reverse=True):
    frequentUNW.append(value[1]) # 部分構造
    frequentUNWNames.append(key)

訳注(04/2020)

Pythonやプログラミング初心者仲間のための蛇足・・・
for ele in highLightFrameUNW.valuesの部分ですが、valuesはDataFrameの実際のデータの値で、NumPyの配列ndarrayです。 DataFramehighLightFrameUNWを参照していただければわかると思いますが、ele[2]は２列目unwantedID列を、ele[3]は3列目unwSubstr列の要素を指定しています。
unwCounts.get(unw, "empty")の部分ですが、get(key[, default])は辞書型においてkey が辞書にあれば key に対する値を、そうでなければ default を返します。ですので、今回は辞書unwCountsにkeyunw(unwntedID, ele[2])があるかどうかを検索し、なければdefalutとしてempytを返すというようになっています。辞書unwCountsのkeyはunwntedID、valueはリスト[出現回数, 好ましくない部分構造]となっているので、辞書に含まれていない好ましくない部分構造の場合は新しくエントリーを辞書に追加して値を出現回数を1に設定、すでに辞書に含まれている構造の場合は出現回数[unw][0]をインクリメントするという手順になっているのだと思います。

訳注ここまで

頻度の高い部分構造の上位８個を描画します。

# 上位８個の頻度の高い部分構造
Draw.MolsToGridImage(mols=list(frequentUNW[0:8]), 
                     legends=list(frequentUNWNames[0:8]),
                     molsPerRow=4)

f:id:magattaca:20200418215914p:plain

ディスカッション

このトークトリアルでは、好ましくない部分構造の検索を行う２つの方法を学びました。: * RDKitにすでに実装されている FilterCatalogクラス、そして、 * RDKitに含まれていない外部のリストとHasSubstructMatch()関数です。

実際のところ、PAINSを部分構造検索で実装することもできます。また、BrenkらによるリストもすでにRDKitに実装されています。実装済みのその他のリストについては (RDKit Documentation)を参照してください。

これまで、HasSubstructMatch()関数を使用してきましたが、これは一つの化合物あたり一つのマッチ結果しか出しません。GetSubstructMatches()　関数を使用すると、一つの化合物に含まれる全ての部分構造を見つけることができます。 同様にPAINSに関して、一つの分子あたり、最初にマッチしたもの GetFirstMatch() だけをみてきました。もし、全てのPAINSをフィルタリングして除去したいのであればこれで十分です。ですが、ある分子のもつ全ての危険な部分構造を見るために、GetMatches()を使用することもできます。一つの分子あたり全てのマッチする部分構造を考慮しているわけではないので、最後に描画した部分構造が実際に最も頻度の高いものであったと言うことはできません。ですが、頻度が高いということだけは明らかです。

見つかった部分構造は２つの異なる方法で処理できます。 * 部分構造検索をフィルターとして適用し、引っかかった化合物を、資金と時間の節約のため、さらなる試験からは除外する * あるいは警告として使用することもできます。好ましくない部分構造を有する分子にフラグを立てることができます。（例えば、化学者や毒性学者•・・といった）専門家の目から見て、経験に基づき判断することができるかもしれません。もし、各部分構造がそこまでクリティカルなものでなければ、スクリーニングの対象化合物として含めることも可能です。

ここでは、機械学習に使用するために、あまりに多くの化合物を削除することはしたくなかったので、好ましくない部分構造をフィルタリングで除去することはしません。また、部分構造によるフィルタリングは、後ほど実際のスクリーニング実験を行うまえに適用することもできます。警告（アラート）のフラッグを設定することも可能で、そうすれば（PAINSやBrenk等どのようなリストにでも順じて）好ましくない部分構造についての情報を保持し、あとで考慮することが可能です。

クイズ

  • なぜスクリーニングライブラリーから「PAINS」を取り除くことを考えるべきなのでしょうか？これらの化合物の問題点とは何でしょうか？
  • ある好ましくない部分構造を取り除く必要がないような状況を見つけることはできますか？
  • このチュートリアルで使用した部分構造はどうやってエンコードされていましたか？

訳以上

追記

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

今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。

PAINSについて:
* 化学の新しいカタチ：RDKitのPAINSフィルターで化合物をスクリーニング

RDKitを用いた部分構造検索
* 化学の新しいカタチ：RDKitを用いた部分構造検索とMCSアルゴリズム

さて、今回の記事では医薬品として好ましくない部分構造を検索し、データセットから取り除く方法が取り上げられています。一方で、好ましくない部分構造が状況によって変わりうるののだ、ということも指摘されています。例えば冒頭で好ましくない部分構造の例としてニトロ基や2-ハロピリジンといったものが指摘されています。ではこれらの構造が「絶対に医薬品としてありえないか？」というとその答えはNOです。

ニトロ基を有する医薬品の例として高血圧・狭心症治療剤のニフェジピン、より最近の承認薬としては大塚製薬の肺結核治療薬デラマニドが挙げられます。

「ニトロ基がなぜ嫌われるのか?」というと生物学者の方にはアッセイの検出を思い浮かべていただくとわかりやすいかもしれません。学生実習の題材でも使われるアルカリホスファターゼの酵素活性の検出では、酵素の機能により放出されるp-nitrophenyl phosphate (pNPP)の蛍光性が検出原理として使われています。光を吸収、呈色するような分子構造、なんとなく飲みたくないですよね。（大学受験を終えたばかりでこんなブログにたどりついてしまった方にはキサントプロテイン反応というと分かりやすいかもしれません。）

また、2-ハロピリジンそのものではありませんが現在話題のファビピラル(アビガン)も類似の部分構造を有しています。ピラジン骨格ですがN横の脱離基（ハロゲン）としては反応性が高い構造かもしれません。「部分構造をどの範囲までとるか？」という基準で議論が分かれそうです。ぜひ専門家の皆様のご意見を伺いたいところです。

TeachOpenCADD トピック2 〜指標による化合物フィルタリング〜

TeachOpenCADD

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

取り上げられている内容の概観のため目次を貼っておきます。

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

github.com

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

以下、訳。

トークトリアル 2

化合物フィルタリング: ADME とリードライクネスのクライテリア

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

Michele Ritschel and Mathias Wajnberg

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

ChEMBLから取得した化合物 (talktorial 1)をリードライクネスのクライテリアに基づきフィルタリングします。我々のスクリーニングライブラリーから薬らしくない分子を取り除くのが目的です。

  • 化合物の生物学的利用可能性（バイオベイラビリティ）に関連するパラメータを計算(リピンスキーのルールオブファイブ)
  • ルールオブファイブのクライテリアに基づきChEMBLから集めた化合物をフィルタリング
  • レーダーチャート形式でパラメーターをプロット

学習の目標

理論

  • ADME - 吸収: absorption, 分布: distribution, 代謝: metabolism and 排出: excretion
  • リードライクネスとリピンスキーのルールオブファイブ
  • リードライクネスの文脈におけるレーダーチャートのバリエーションと解釈

実践

  • 例となる化合物について物理化学的パラメータを計算
  • 複数の分子についてそれぞれの物理化学的パラメーターを比較するため棒グラフを作成
  • ルールオブファイブを満たしているかチェックするための関数を定義
  • ChEMBLから取得したデータセット全体にルールオブファイブを適用
  • ルールオブファイブでフィルターをかけたデータセットのレーダーチャートを作成。これにより一つのプロットでルールオブファイブのクライテリアに関するプロパティの可視化が可能

レファレンス

理論

バーチャルスクリーニングでは、ある化合物が特定のターゲット分子に対して結合し、相互作用するかということを予測することができます。ですが、新しい医薬品を見つけ出すためには、その化合物がターゲット分子に到達し、最終的には体内から好ましい経路で取り除かれることが重要です。したがって、化合物が実際に体の中に取り込まれ、ターゲット分子に到達するためのて特定のバリアーを通過することができるかということも考慮しなければいけません。代謝的に安定か？そしてターゲット分子に対して作用しなくなったらどのように排出されるのか？これらの過程は薬物動態学（pharmacokinetics）の分野で研究されています。薬物動力学（pharmacodynamics）（「医薬品は私たちの体に何をするのか？」）とは対照的に、薬物動態学は「私たちの体のなかで医薬品に対して何が起こるのか？」という問いを扱います。

ADME

薬物動態学は主に４つのステップに分けられます。 Absorption：吸収、 Distribution：分布、 Metabolism：代謝、そして Excretion：排出です。 これら４ステップは ADME　と略されます。ADME(T) はToxicology（毒性学）も含むことがあります。 以下ではADMEの過程がより詳細に議論されています。 (ADME wikipedia and Mol Pharm. (2010), 7(5), 1388-1405)

  • 吸収：ある物質が体内に取り込まれる量と取り込みに必要な時間は、各個人とその状態によって異なる様々な要因と、その物質の特性に依存しています。（低い）化合物溶解度と胃内容排出時間（gastric emptying time）、腸通過時間（intestinal transit time）、胃の中での化合物（不）安定性、そして腸壁の（不）透過能、といった要因が全て、医薬品が、例えば経口投与、吸入、そして皮膚への接触のあとでどの程度吸収されるかに影響をあたる可能性があります。

  • 分布：吸収された物質の分布、即ち体内における、血液とそれぞれの組織の間の分布、血液脳関門の透過は、局所的な血流速度や化合物のサイズと極性、そして血清中のタンパク質やトランスポーター酵素への結合の影響を受けます。毒性学における重要な要因としては脂肪組織への疎水性の高い物質の蓄積、あるいは血液脳関門の透過があげられます。

  • 代謝：化合物が体内に入るとすぐに、通常、化合物の代謝が始まります。つまり実際には、この化合物の一部だけがターゲット分子に到達します。主に肝臓と腎臓の酵素が生体異物（xenobiotics：体の外からきた物質）の分解を担っています。毒性のある化合物が取り除かれる場合、吸収された物質の量を減らすことは好ましいことです。その一方で、化学物質の変換は新たな毒性代謝物をうみだす可能性もあります。

  • 排出：化合物と代謝物は排出によって体内から取りのぞかれる必要があります。通常、腎臓（尿）あるいは糞によって排出されます。不完全な排出は結果として、外来物質の蓄積、あるいは通常の代謝を妨げることに繋がります。
Figure 1:ヒト体内におけるADME過程
(figure taken from openclipart.org and adapted)

リードライクネスとリピンスキーのルールオブファイブ

リード化合物 とは有望な特性をもった医薬品開発の候補となる化合物です。望ましい医薬品を見つけ出すために構造修飾を行う出発点となる構造です。生理活性（目標とするターゲット分子に結合する化合物）に加えて、好ましいADME特性が有効な医薬品をデザインする上で重要なクライテリアです。

生物学的利用可能性（バイオアベイラビリティ）は重要なADME特性で、化合物の構造のみに基づいてこの特性を計るために、リピンスキーのルールオブファイブ（Lipinski's rule of five）が発明されました。これは経験則で、化合物の経口投与によるバイオアベイラビリティを見積もるのに役立ちます。

ルールオブファイブによると、ある化合物が次のルールの一つ以上を破っているとき、経口で体内に吸収される可能性が低くなるとされいます。そのルールとは即ち：

  • 分子量が500ダルトン以下
  • 水素結合アクセプターが10以下
  • 水素結合ドナーが５以下
  • LogP (オクタノールー水　分配係数) <= 5

LogP は分配係数あるいはオクタノールー水分配係数とも呼ばれます。通常、疎水性（例：1-オクタノール）と親水性（例：水）の相の間における化合物の分配を測定します。

疎水性分子は水に対する溶解度が低くなる一方で、より親水性の分子（例：水素結合のアクセプターとドナーの数が多いもの）や大きな分子（分子量の大きいもの）は、リン脂質二重膜を透過しづらくなります。

ルールオブファイブに関して、全ての数字が５の倍数になっています。これがこのルールの名称の起源となっています。

(Adv. Drug Deliv. Rev. (1997), 23, 3-25)

レーダーチャート

ルールオブファイブに関する分子の特性を計算した後で、その特性を可視化できると役に立ちます。Ritchieら(Drug. Discov. Today (2011), 16(1-2), 65-72)　はADMEに関する特性のグラフィカルな表現についての概観を示しています。分子の特性を可視化し、創薬化学者による解釈を手助けするために様々な手法があります（例：Craigプロット、フラワープロット、黄金の三角地帯（golden triangle））。

このチュートリアルではPythonの作図ライブラリである matplotlib を使ってレーダーチャートを作成する方法を学びます。その見た目から、レーダーチャート(radar charts wikipedia) はしばしばスパイダープロット、あるいはクモの巣グラフと呼ばれます。３６０度に沿って配置され、それぞれの条件に対して、中心から出発する軸を一つ持っています。各パラメーターの値は軸に沿ってプロットされ直線で結ばれています。影のついたエリアはパラメーターが条件を満たす領域を示します。

Figure 2: 化合物データセットの物理化学的特性を表すレーダーチャート

実践

化合物例を定義し可視化する

ChEMBLから取り出したデータセット全体で実践する前に、４つの化合物を例として取り上げ、その化学的特性を調査したいと思います。必要なライブラリをインポートし、４化合物のSMILESを使って描画します。

from rdkit import Chem
from rdkit.Chem import Descriptors
import pandas as pd
from rdkit.Chem import Draw
import numpy as np

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.lines import Line2D
from math import pi

smiles_1 = 'CCC1C(=O)N(CC(=O)N(C(C(=O)NC(C(=O)N(C(C(=O)NC(C(=O)NC(C(=O)N(C(C(=O)N(C(C(=O)N(C(C(=O)N(C(C(=O)N1)C(C(C)CC=CC)O)C)C(C)C)C)CC(C)C)C)CC(C)C)C)C)C)CC(C)C)C)C(C)C)CC(C)C)C)C' # Cyclosporine
smiles_2 = 'CN1CCN(CC1)C2=C3C=CC=CC3=NC4=C(N2)C=C(C=C4)C' # Clozapine
smiles_3 = 'CC1=C(C(CCC1)(C)C)C=CC(=CC=CC(=CC=CC=C(C)C=CC=C(C)C=CC2=C(CCCC2(C)C)C)C)C' # Beta-carotene
smiles_4 = 'CCCCCC1=CC(=C(C(=C1)O)C2C=C(CCC2C(=C)C)C)O' # Cannabidiol
smiles_list = [smiles_1, smiles_2, smiles_3, smiles_4]
names_list = ['cyclosporine', 'clozapine', 'beta-carotene', 'cannabidiol']
mol_list = [Chem.MolFromSmiles(smiles) for smiles in smiles_list]

Draw.MolsToGridImage(mol_list, legends=names_list, molsPerRow=4)

f:id:magattaca:20200412172944p:plain

ルールオブファイブに関する分子の特性を計算しプロットする

ルールオブファイブに関する化学特性を計算し視覚的に比較します：

MWs = [Descriptors.ExactMolWt(mol) for mol in mol_list]
HBAs = [Descriptors.NumHAcceptors(mol) for mol in mol_list]
HBDs = [Descriptors.NumHDonors(mol) for mol in mol_list]
LogPs = [Descriptors.MolLogP(mol) for mol in mol_list]
parameters = [MWs, HBAs, HBDs, LogPs]
print('Molecular weight of the four compounds:',MWs)
#  Molecular weight of the four compounds: [1201.841367992, 306.184446704, 536.438201792, 314.2245802]
  • 分子ごとの特性を棒グラフでプロット
# 2x2のプロットの枠を作成
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2)
axes = [ax1, ax2, ax3, ax4]
x = np.arange(1, len(mol_list)+1)
colors = ['red', 'green', 'blue', 'cyan']

# サブプロットを作成
for index in x-1:
    axes[index].bar(x, parameters[index], color=colors)

# ルールオブファイブの閾値を点線で追加
ax1.axhline(y=500, color="black", linestyle="dashed")
ax1.set_title("molecular weight (Da)")
ax2.axhline(y=10, color="black", linestyle="dashed")
ax2.set_title("# h-bond acceptors")
ax3.axhline(y=5, color="black", linestyle="dashed")
ax3.set_title("# h-bond donors")
ax4.axhline(y=5, color="black", linestyle="dashed")
ax4.set_title("logP")

# 凡例の追加
legend_elements = [mpatches.Patch(color=colors[i], label=names_list[i]) for i in range(len(mol_list))]
legend_elements.append(Line2D([0], [0], color="black", ls="dashed", label="Threshold"))
fig.legend(handles=legend_elements, bbox_to_anchor=(1.25, 0.5))

# サブプロットと凡例を図に追加
plt.tight_layout()

plt.show()

f:id:magattaca:20200412173052p:plain

上の棒グラフでは４例の化合物についてルールオブファイブの特性 （分子量、水素結合ドナー数、水素結合アクセプター数、LogP）を比較しました。４例の医薬品分子が異なる特性を持つことが見てとれます。次のステップでは各化合物のそれぞれについてルールオブファイブを逸脱しているかどうか調べます。

リピンスキーのルールオブファイブに従っているかの調査

ある化合物がルールオブファイブを逸脱しているか否かを調べるための関数を定義し、４つの化合物例に適用します。

def test_rule_of_five(smi):
    m = Chem.MolFromSmiles(smi)
    
    # ルールオブファイブの化学特性を計算
    MW = Descriptors.ExactMolWt(m)
    HBA = Descriptors.NumHAcceptors(m)
    HBD = Descriptors.NumHDonors(m)
    LogP = Descriptors.MolLogP(m)
    
    # ルールオブファイブの条件
    conditions = [MW <= 500, HBA <= 10, HBD <= 5, LogP <= 5]
    # ４つの条件のうち逸脱が一つ以下の時にTrueを返す
    return conditions.count(True) >= 3

for i in range(len(smiles_list)):
    smi=smiles_list[i]
    name=names_list[i]
    print("Rule of five accepted for %s: %s "%(name,test_rule_of_five(smi)))                

# Rule of five accepted for cyclosporine: False 
#  Rule of five accepted for clozapine: True 
#  Rule of five accepted for beta-carotene: False 
#  Rule of five accepted for cannabidiol: True

我々の test_rule_of_five 関数によると４つの化合物例のうち２つはルールオブファイブを満たしていません。ここから、シクロスポリン（cyclosporin）とβ-カロテン（betacarotene）は経口で体内に吸収される可能性が低いと解釈できます。４つとも全て医薬品として市場で手に入るので、それぞれのターゲット分子にどうにかして到達しているはずです。ルールの例外である可能性もありますし、経口投与とは異なる投与ルートで使われている可能性もあります。

EGFRデータセットにルールオブファイブを適用

test_rule_of_five 関数を使って、リピンスキーのルールオブファイブに従っているかを基準にして、目的のメインデータセットにフィルターをかけます。

  • ルールオブファイブに関する化学特性を全て返すように関数を調整
  • メインとなるデータフレームを読み込み (ChEMBL_df)
  • ルールオブファイブの関数をChEMBL_df　に適用
  • ルールを２つ以上逸脱しているかによって ChEMBL_df にフィルターをかける
  • フィルタリングした後のデータフレームを保存
def df_rule_of_five(df):
    
    smi = df['smiles']
    m = Chem.MolFromSmiles(smi)
    
    # ルールオブファイブの化学特性を計算
    MW = Descriptors.ExactMolWt(m)
    HBA = Descriptors.NumHAcceptors(m)
    HBD = Descriptors.NumHDonors(m)
    LogP = Descriptors.MolLogP(m)
    
    # ルールオブファイブの条件
    conditions = [MW <= 500, HBA <= 10, HBD <= 5, LogP <= 5]
    
    # 特性の値とルールオブファイブに従っているか否かについての情報をもつpandasの行を作成
    return pd.Series([MW, HBA, HBD, LogP, 'yes']) if conditions.count(True) >= 3 else pd.Series([MW, HBA, HBD, LogP, 'no'])

ChEMBL_df = pd.read_csv('../data/T1/EGFR_compounds.csv', index_col=0)
print(ChEMBL_df.shape)
ChEMBL_df.head()
# (5428, 5)
molecule_chembl_id units IC50 smiles pIC50
0 CHEMBL63786 nM 0.003 Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 11.522879
1 CHEMBL53711 nM 0.006 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849
2 CHEMBL35820 nM 0.006 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849
3 CHEMBL53753 nM 0.008 CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.096910
4 CHEMBL66031 nM 0.008 Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 11.096910 
# ルールオブファイブの結果を得るために関数を適用 (しばらく時間がかかる可能性があります)
rule5_prop_df = ChEMBL_df.apply(df_rule_of_five, axis=1)

# 条件列に名前をつける
rule5_prop_df.columns= ['MW', 'HBA', 'HBD', 'LogP', 'rule_of_five_conform']

# 計算した値をもつデータセットを連結（concatenate）
ChEMBL_df = ChEMBL_df.join(rule5_prop_df)

# 空の行を削除 --> rule of five
filtered_df = ChEMBL_df[ChEMBL_df['rule_of_five_conform']=='yes']

# データに関する情報
print('# of compounds in unfiltered data set:', len(ChEMBL_df))
print('# of compounds in filtered data set:', len(filtered_df))
print("# of compounds not compliant with Lipinski's rule of five:", (len(ChEMBL_df)-len(filtered_df)))

# フィルターをかけたデータの保存
filtered_df.to_csv('../data/T2/EGFR_compounds_lipinski.csv', sep=';') 
filtered_df.head(10)

#  # of compounds in unfiltered data set: 5428
#  # of compounds in filtered data set: 4523
#  # of compounds not compliant with Lipinski's rule of five: 905
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.28910 yes
1 CHEMBL53711 nM 0.006 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849 343.043258 5 1 3.59690 yes
2 CHEMBL35820 nM 0.006 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849 387.058239 5 1 4.93330 yes
3 CHEMBL53753 nM 0.008 CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.096910 329.027607 5 2 3.57260 yes
4 CHEMBL66031 nM 0.008 Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 11.096910 339.011957 4 2 4.01220 yes
5 CHEMBL176582 nM 0.010 Cn1cnc2cc3ncnc(Nc4cccc(Br)c4)c3cc21 11.000000 353.027607 5 1 4.02260 yes
6 CHEMBL174426 nM 0.025 Cn1cnc2cc3c(Nc4cccc(Br)c4)ncnc3cc21 10.602060 353.027607 5 1 4.02260 yes
7 CHEMBL29197 nM 0.025 COc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OC 10.602060 359.026939 5 1 4.15310 yes
8 CHEMBL1243316 nM 0.030 C#CCNC/C=C/C(=O)Nc1cc2c(Nc3ccc(F)c(Cl)c3)c(C#N... 10.522879 477.136781 6 3 4.75878 yes
9 CHEMBL3613702 nM 0.037 C=CC(=O)Nc1ccc2ncnc(Nc3cc(F)c(Cl)c(Cl)c3)c2c1 10.431798 376.029395 4 2 4.94380 yes

ルールオブファイブの特性を可視化するためのレーダープロット

まずデータセットの平均と標準偏差を求める関数を定義します。

これらの統計量はあとでデータセットのリピンスキーのルールオブファイブに関するパラメーターをプロットすることに使います。

def get_properties_stats(data_df):
    """
    データセットの物理化学特性の平均と標準偏差を計算する関数
    
    Input: 
    化合物ごとの物理化学特性の値、 HBD, HBA, MW そして　LogPを列として
    （正確にこのままの列名で）含むデータセット
    
    Output:
    物理化学特性（行）それぞれについて平均と標準偏差の値（列）をもつデータフレーム
    """
    properties = ["HBD", "HBA", "MW", "LogP"]
    
    data_stats = []
    
    for i in properties:
        std = data_df[i].std()
        mean = data_df[i].mean()
        da = pd.DataFrame([[round(mean, 2), round(std, 2)]], index=[i], columns=["mean", "std"])
        data_stats.append(da)
    
    data_stats = pd.concat(data_stats)
    
    return data_stats

リピンスキーのルールオブファイブに一致する化合物データセットの統計量を計算します（フィルタリングしたデータセット）。

stats_rof = get_properties_stats(filtered_df)
stats_rof
mean std
HBD 1.88 1.01
HBA 5.97 1.88
MW 413.03 87.93
LogP 4.07 1.19

リピンスキーのルールオブファイブに合致しない化合物データセットの統計量を計算します。

stats_not_rof = get_properties_stats(ChEMBL_df[ChEMBL_df['rule_of_five_conform']=='no'])
stats_not_rof
mean std
HBD 2.30 1.74
HBA 7.99 2.37
MW 588.96 103.05
LogP 5.96 1.44

レーダーチャートで化合物特性を可視化する関数を作成します。stackoverflowのチュートリアルに従って行います。

def plot_radarplot(data_stats, output_path):
    """
    ４つの物理化学特性（HBD、HBA、MW、LogP)の平均と標準偏差に基づくレーダープロットを作成する関数
    
    Input: 
    それぞれの物理化学特性（行）について平均と標準偏差（列）を有するデータフレーム。
    
    Output:
    レーダープロット (ファイルとして保存し、Jupyter notebookにも表示する)。
    """

    # 線を引くためのデータポイントを取得
    std_1 = [data_stats["mean"]["HBD"] + data_stats["std"]["HBD"], 
             (data_stats["mean"]["HBA"]/2) + (data_stats["std"]["HBA"]/2), 
             (data_stats["mean"]["MW"]/100) + (data_stats["std"]["MW"]/100), 
             data_stats["mean"]["LogP"] + data_stats["std"]["LogP"]]
    std_2 = [data_stats["mean"]["HBD"] - data_stats["std"]["HBD"], 
             (data_stats["mean"]["HBA"]/2) - (data_stats["std"]["HBA"]/2), 
             (data_stats["mean"]["MW"]/100) - (data_stats["std"]["MW"]/100), 
             data_stats["mean"]["LogP"] - data_stats["std"]["LogP"]]
    mean_val = [data_stats["mean"]["HBD"], (data_stats["mean"]["HBA"]/2), 
                (data_stats["mean"]["MW"]/100), data_stats["mean"]["LogP"]]

    # （塗りつぶされた）領域のデータポイント（ルールオブファイブ）を取得
    rule_conditions = [5, (10/2), (500/100), 5]
    
    # 特性の名称を定義
    parameters = ['# H-bond donors', '# H-bond acceptors/2', 'Molecular weight (Da)/100', 'LogP']

    # 
    N = len(rule_conditions)

    # フォントサイズを指定
    fontsize = 16

    # 条件の軸の角度
    x_as = [n / float(N) * 2 * pi for n in range(N)]

    # 円型チャートなので各リストの最初の値のコピーを、
    # 各リストの最後にデータとともに追加する必要があります。
    std_1 += std_1[:1]
    std_2 += std_2[:1]
    mean_val += mean_val[:1]
    rule_conditions += rule_conditions[:1]
    x_as += x_as[:1]

    # 図のサイズを指定
    plt.figure(figsize=(8,8))

    # 軸の色を指定
    plt.rc('axes', linewidth=2, edgecolor="#888888")

    # 極座標プロット（polar plot）の作成
    ax = plt.subplot(111, polar=True)

    # 時計回りを指定。つまり：
    ax.set_theta_offset(pi / 2)
    ax.set_theta_direction(-1)

    # y軸のラベル位置を指定
    ax.set_rlabel_position(0)

    # グリッドの色と線のスタイルを指定
    ax.xaxis.grid(True, color="#888888", linestyle='solid', linewidth=2)
    ax.yaxis.grid(True, color="#888888", linestyle='solid', linewidth=2)

    # 放射状の数の指定とラベルの除去
    plt.xticks(x_as[:-1], [])

    # y軸の印（ytick）の指定
    plt.yticks([1, 3, 5, 7], ["1", "3", "5"], size=fontsize)

    # 軸の範囲（limit）の指定
    plt.ylim(0, 7)

    # データをプロット
    # 平均値
    ax.plot(x_as, mean_val, 'b', linewidth=3, linestyle='solid', zorder=3)

    # 標準偏差
    ax.plot(x_as, std_1, 'm', linewidth=2, linestyle='dashed', zorder=3, color='#111111')
    ax.plot(x_as, std_2, 'y', linewidth=2, linestyle='dashed', zorder=3, color='#333333')

    # 領域を塗りつぶす
    ax.fill(x_as, rule_conditions, "#3465a4", alpha=0.2)

    # 正しくフィットしていることを確かめるためにy軸の印のラベルを描く
    for i in range(N):
        angle_rad = i / float(N) * 2 * pi
        if angle_rad == 0:
            ha, distance_ax = "center", 1
        elif 0 < angle_rad < pi:
            ha, distance_ax = "left", 1
        elif angle_rad == pi:
            ha, distance_ax = "center", 1
        else:
            ha, distance_ax = "right", 1
        ax.text(angle_rad, 7 + distance_ax, parameters[i], size=fontsize,
                horizontalalignment=ha, verticalalignment="center")

    # 左上のプロットに合わせて凡例を追加 
        labels = ('Mean', 'Mean + std', 'Mean - std', 'Rule of five area')
        legend = ax.legend(labels, loc=(1.1, .7),
                           labelspacing=0.3, fontsize=fontsize)
    plt.tight_layout()

    # プロットを保存　- テキストボックスを含めるために bbox_inches　を使用：
    # https://stackoverflow.com/questions/44642082/text-or-legend-cut-from-matplotlib-figure-on-savefig?rq=1
    plt.savefig(output_path, dpi=300, bbox_inches="tight", transparent=True)

    # 極座標プロットを表示
    plt.show()

まず、ルールオブファイブでフィルタリングしたデータセットをプロットします。

plot_radarplot(stats_rof, "../data/T2/radarplot_rof.png")

f:id:magattaca:20200412173555p:plain

上で作られたレーダチャートで、青色の四角は物理化学特性がルールオブファイブの基準の中に収まっている領域を表します。青線はフィルタリングしたデータセットの平均値を結んだもので、点線は標準偏差を結んだものです。平均値はリピンスキーのルールのどれも逸脱していないということがわかります。ですが、標準偏差をみると、いくつかの特性値はまだ閾値よりも大きくなっています。これは許容範囲内です。４つの特性のルールのうち一つは逸脱する可能性があるということを覚えておいてください。

次に、ルールオブファイブを逸脱した化合物を眺めます。

plot_radarplot(stats_not_rof, "../data/T2/radarplot_not_rof.png")

f:id:magattaca:20200412173624p:plain

ルールオブファイブに従っていない化合物の特性のほとんどが、LogPと分子量によってルールを外れていることがわかります。

ディスカッション

リピンスキーのルールオブファイブは経口吸収性（バイオアベイラビリティ）に焦点をあてています。医薬品は経口以外の投与ルートでも投与可能です（例：吸入、経皮、注射）。ルールオブファイブは経口投与における生物学的利用可能性の指標であって、例外はあるのだということに注意してください。バイオアベイラビリティを題材に、いくつかのADME特性の一つを見てきました。

ADME特性の全体像を手にするために利用可能なwebサーバーとプログラムがあります。　例： SwissADME.

クイズ

  • ルールオブファイブで表現された化学特性はADMEに対してどのような影響を与えうるでしょうか？
  • ３つ、あるいは４つのルールを逸脱するような分子を見つけるか、デザインしてください。
  • 上で作成したレーダーチャートに新しい化合物の情報を追加でプロットするにはどうすればよいでしょうか？

訳以上

追記

今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
ところでレーダーチャートの利用と書き方を初めて見ました。
設定項目が多くて面倒ですが一つのプロットで複数のパラメーターが確認できるという点では便利そうですね。

今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。

Drug likenessの指標について:
* 化学の新しいカタチ：色々な薬らしさの指標
* Druglikenessについてのよもやま話

RDKitによる記述子計算:
* 化学の新しいカタチ：RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ

TeachOpenCADD トピック1 〜化合物データの取得 (ChEMBL)~

TeachOpenCADD

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

取り上げられている内容の概観のため目次を貼っておきます。

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

github.com

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

以下、訳。

トークトリアル１

化合物データの取得 (ChEMBL)

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

Paula Junge and Svetlana Leng

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

ChEMBLからデータを抽出する方法の学習:

  • ある特定の標的に対して評価済みのリガンドを見つける
  • 取得可能な生理活性データでフィルタリング
  • pIC50値の計算
  • データフレームを結合し、取り出した分子を描画

学習の目標

理論

実践

目標：与えられた標的に対する生理活性データのある化合物のリストを取得すること

  • ChEMBLデータベースへの接続
  • 標的（ターゲット分子）に関するデータの取得 (EGFR kinase)
  • 生理活性データ
    • 生理活性をダウンロードしフィルタリングする
    • データの整型と変換
  • 化合物データ
    • 化合物リストの取得
    • 出力データの準備
  • 出力
    • 最もpIC50値の大きい分子を描画
    • アウトプットファイルの書き出し

レファレンス

理論

ChEMBLデータベース

  • 大規模な生理活性の公開データベース
  • 現在のデータ内容 (10.2018時点):
    • 1.8 百万以上のユニークな構造の化合物
    • 1 百万アッセイに基づく 15 百万以上の活性値
    • ∼12 000 のターゲット分子に紐づけられたアッセイ
  • データソース に含まれているもの：科学文献、PubChemバイオアッセイ、顧みられない病気のための創薬イニシアティブ(Drugs for Neglected Diseases Initiative (DNDi))、BindingDB データベース、などなど
  • ChEMBLのデータは Webインターフェースあるいは EBI-RDFプラットフォーム 、 ChEMBL Webサービス　からアクセスすることができます。

ChEMBL Webサービス

f:id:magattaca:20200411143208j:plain

Figure 1: "ChEMBL webサービス スキーマダイアグラム。楕円型はChEMBL webサービスのリソースを示し、二つのリソースを繋ぐラインは、それらが共通の特徴を共有していることを示します。矢印の向きはリソースタイプについての一次情報を見つけられる場所を示しています。点線は２つのリソースの関係が異なる形式で振舞うことを示しています。例えば、Image リソースは Molecule のグラフィカルな描写を提供します。Figureと説明は次の文献からとりました: Nucleic Acids Res. (2015), 43, 612-620.

ChEMBL webリソースクライアント

  • ChEMBLのデータにアクセスするためのPythonクライアントライブラリ
  • HTTPSプロトコルとのやり取りを行う
  • 結果の遅延評価（lazy evaluation） -> ネットワークリクエストを減らす

化合物活性評価

IC50

Figure 2: IC50値を算出する方法を視覚的に表した図:　阻害を縦軸、log(濃度)を横軸にデータを配置し、阻害の最大と最小を決定する。曲線が50%阻害レベルと交差する点の濃度がIC50

pIC50

  • IC50値の比較を容易にするため、pIC50値をlogスケールで次のように定義しています
    $pIC{50} = -log{10}(IC{50}) $ （$ IC{50}$の単位はM）
  • pIC50が大きくなるほど指数関数的に薬効が強くなっていることを示します
  • pIC50はモル濃度の指標で与えられます (mol/L or M)
    • IC50はM単位に換算したあとで、pIC50に変換する必要があります。
    • nMの場合は: $pIC{50} = -log{10}(IC{50}*10^{-9})= 9-log{10}(IC_{50}) $

IC50と、pIC50に加えて、その他の生理活性の指標として、平衡定数 KI と50%効果濃度(half maximal effective concentration (EC50 ))といったものが使われます。

実践

それでは、我々の興味のあるターゲット分子であるEGFRキナーゼに対して、評価が実施された全ての分子をダウンロードしたいと思います。

ChEMBLデータベースへの接続

まず、ChEMBL webリソースクライアントと他のpythonライブラリーをインポートします。

from chembl_webresource_client.new_client import new_client
import pandas as pd
import math
from rdkit.Chem import PandasTools

訳注(04.2020時点) from chembl_webresource_client.new_client import new_client を実行すると以下のようなでエラーが出ました。

Exception: svg+xml is not an available format (['xml', 'json', 'yaml', 'svg'])

こちらの問題はclientのバージョンを上げることで解決できます。イシュー #71を参照してください。アップデートは

pip install -U chembl-webresource-client

で実行可能でした。 「-U」を忘れずに！
訳注ここまで

APIアクセスのためのリソースオブジェクトを作成します

targets = new_client.target
compounds = new_client.molecule
bioactivities = new_client.activity

ターゲット分子のデータ

  • UniProt Webサイト (https://www.uniprot.org/)　から目的とするターゲット分子（EGFRキナーゼ）の UniProt-ID (http://www.uniprot.org/uniprot/P00533) を取得します。
  • UniProt-ID を使用してターゲット分子の情報を取得します。
  • 他のターゲット分子に興味があるときは、ターゲットに合わせて異なるUniProt-IDを選択してください
uniprot_id = 'P00533'
# ChEMBLからターゲット分子の情報を取得しますが、
# 取得する項目を特定のものに絞って取得します。
target_P00533 = targets.get(target_components__accession=uniprot_id) \
                       .only('target_chembl_id', 'organism', 'pref_name', 'target_type')
print(type(target_P00533))
pd.DataFrame.from_records(target_P00533)

# <class 'chembl_webresource_client.query_set.QuerySet'>
cross_references organism pref_name species_group_flag target_chembl_id target_components target_type tax_id
0 [{'xref_id': 'O43451', 'xref_name': None, 'xre... Homo sapiens Maltase-glucoamylase False CHEMBL2074 [{'accession': 'O43451', 'component_descriptio... SINGLE PROTEIN 9606.0
1 NaN Homo sapiens Epidermal growth factor receptor erbB1 NaN CHEMBL203 NaN SINGLE PROTEIN NaN
2 NaN Homo sapiens Epidermal growth factor receptor and ErbB2 (HE... NaN CHEMBL2111431 NaN PROTEIN FAMILY NaN
3 NaN Homo sapiens Epidermal growth factor receptor NaN CHEMBL2363049 NaN PROTEIN FAMILY NaN
4 NaN Homo sapiens MER intracellular domain/EGFR extracellular do... NaN CHEMBL3137284 NaN CHIMERIC PROTEIN NaN

訳注(0４.2020時点)
出力したtypeはQuerySetとなっています。Django QuerySetインターフェースに基づいたデザインとなっているそうです。詳しくはchemble_webresource_clientのGitHubを参照してください 出力結果はオリジナルのJupyter Notebookと異なっています。

  • only(~) で項目を絞っているはずですが、データフレームに表示された項目は絞られていません。
  • また、chembl_webresouce_clientのonlyオペレーターの使用方法では引数をリストにすることが推奨されています。
  • index 0 のエントリーが「CHEMBL2074」 となっています。こちらはEGFRではなさそうです。エラーでしょうか？

参考までに、データの取得と絞り込みを２段階に分けた方法を下に示します。データフレームのカラムを選択しているだけです。

target_P00533 = targets.get(target_components__accession=uniprot_id) 
selected_target_values = ['organism', 'pref_name', 'target_chembl_id', 'target_type']
df_trarget_tmp = pd.DataFrame.from_records(target_P00533)[selected_target_values]
df_trarget_tmp
organism pref_name target_chembl_id target_type
0 Homo sapiens Maltase-glucoamylase CHEMBL2074 SINGLE PROTEIN
1 Homo sapiens Epidermal growth factor receptor erbB1 CHEMBL203 SINGLE PROTEIN
2 Homo sapiens Epidermal growth factor receptor and ErbB2 (HE... CHEMBL2111431 PROTEIN FAMILY
3 Homo sapiens Epidermal growth factor receptor CHEMBL2363049 PROTEIN FAMILY
4 Homo sapiens MER intracellular domain/EGFR extracellular do... CHEMBL3137284 CHIMERIC PROTEIN

訳注ここまで

ヒットしたエントリーを確認したのち、最初のエントリーを目的のターゲット分子として選択

CHEMBL203: 単一のタンパク質（single protein）で human Epidermal growth factor receptor (EGFR, 別名 erbB1) を表しています。

target = target_P00533[0]
target
#    {'organism': 'Homo sapiens',
#    'pref_name': 'Epidermal growth factor receptor erbB1',
#    'target_chembl_id': 'CHEMBL203',
#    'target_type': 'SINGLE PROTEIN'}

訳注(04.2020時点)

DataFrameのindex 0 と同一の結果（「CHEMBL2074の出力」）となるかと思いましたが、オリジナルのGitHub のノートブック同様に「CHEMBL203」のが出力されました。QuerySetをDataFrameに変換するタイミングでエラーがあるのでしょうか？

訳注ここまで

選択したChEMBL-IDを保存します。

chembl_id = target['target_chembl_id']
chembl_id
#  'CHEMBL203'

生理活性データ

それでは、目的とするターゲット分子の生理活性データをクエリとして投げたいと思います。

ターゲット分子の生理活性をダウンロードしてフィルタリング

この工程では、生理活性データをダウンロードし、フィルターをかけます。次の項目のみ考慮します。

  • ヒトタンパク質（human proteins）
  • 生理活性タイプ（bioactivity type IC50）
  • 精確な測定結果（関係式が '='　となっているもの)
  • 結合データ (アッセイタイプ 'B')
bioact = bioactivities.filter(target_chembl_id = chembl_id) \
                      .filter(type = 'IC50') \
                      .filter(relation = '=') \
                      .filter(assay_type = 'B') \
                      .only('activity_id','assay_chembl_id', 'assay_description', 'assay_type', \
                            'molecule_chembl_id', 'type', 'units', 'relation', 'value', \
                            'target_chembl_id', 'target_organism')
# len(bioact), len(bioact[0]), type(bioact), type(bioact[0])
len(bioact[0]), type(bioact), type(bioact[0])

#  (43, chembl_webresource_client.query_set.QuerySet, dict)

訳注(04.2020時点)

len(bioact) はNonTypeのため表示できないといったエラーが出ました。コメントアウトし、残りの３つを表示させています。
len(bioact) に相当するデータ（ヒットした化合物の数）は以降のDataFrame化、サイズ確認の段階で確認できます。

訳注ここまで

ChEMBLデータベースにクエリを投げるのが難しければ、先のセルのクエリの結果（11 April 2019）を含むファイルを提供しています。Pythonのパッケージ pickleを使いました。pickelはPythonのオブジェクトをシリアライズし、ファイルに保存、あとでプログラムに読み込めるようにします。（オブジェクトのシリアライズについて、さらに学びたいときには次を参照してください DataCamp )

pickle化した化合物を読み込むには次のセルのコメントアウトを外して実行してください。

#import pickle
#bioact = pickle.load(open("../data/T1/EGFR_compounds_from_chembl_query_20190411.p", "rb"))

生理活性データを整型し変換する

データは辞書のリストととして保存されています。

bioact[0]

{'activity_comment': None,
 'activity_id': 31863,
 'activity_properties': [],
 'assay_chembl_id': 'CHEMBL663853',
 'assay_description': 'Inhibitory concentration against human DNA topoisomerase II, alpha mediated relaxation of pBR322; no measurable activity',
 'assay_type': 'B',
 'bao_endpoint': 'BAO_0000190',
 'bao_format': 'BAO_0000357',
 'bao_label': 'single protein format',
 'canonical_smiles': 'c1ccc(-c2nc3c(-c4nc5ccccc5o4)cccc3o2)cc1',
 'data_validity_comment': None,
 'data_validity_description': None,
 'document_chembl_id': 'CHEMBL1137930',
 'document_journal': 'Bioorg. Med. Chem. Lett.',
 'document_year': 2004,
 'ligand_efficiency': None,
 'molecule_chembl_id': 'CHEMBL113081',
 'molecule_pref_name': None,
 'parent_molecule_chembl_id': 'CHEMBL113081',
 'pchembl_value': None,
 'potential_duplicate': False,
 'qudt_units': 'http://www.openphacts.org/units/Nanomolar',
 'record_id': 206172,
 'relation': '>',
 'src_id': 1,
 'standard_flag': True,
 'standard_relation': '>',
 'standard_text_value': None,
 'standard_type': 'IC50',
 'standard_units': 'nM',
 'standard_upper_value': None,
 'standard_value': '100000.0',
 'target_chembl_id': 'CHEMBL1806',
 'target_organism': 'Homo sapiens',
 'target_pref_name': 'DNA topoisomerase II alpha',
 'target_tax_id': '9606',
 'text_value': None,
 'toid': None,
 'type': 'IC50',
 'units': 'uM',
 'uo_units': 'UO_0000065',
 'upper_value': None,
 'value': '100.0'}

andasのDataFrameに変換します。（数分かかるかもしれません）

訳注(04.2020時点)
上の出力でもクエリでのデータ取得の段階での項目の絞り込みがうまく機能していません。
オリジナルのノートブックでは次にpandas.DataFrame.from_recordsによりDataFrameに変換していますが、項目数が多すぎるためかうまく変換できませんでした。
そこで冗長ですが、以下のセルでは一度csvファイルに出力し、csvをDataFrameに読み込むという作業を行なっています。
csv出力の段階で項目の絞り込みを実施ししています。（クエリのonly の代替）。
また、DataFrame変換後、項目の値の条件によるフィルタリングを念のため再度実施しています。
訳注ここまで

import csv

# 出力したいフィールド （csvのheaderにする）
fields = ['activity_id','assay_chembl_id', 'assay_description', 'assay_type', 'molecule_chembl_id', 'type', 'units', 'relation', 'value', 'target_chembl_id', 'target_organism']

# bioact.csvというファイルを作成し出力する
with open('../data/T1/bioact.csv', 'w') as f:
    writer = csv.writer(f)
    
    # ヘッダーの書き込み
    writer.writerow(fields)
    
    # activityのデータを１行ずつ書き込み
    for act_data in bioact:
        row = []
        
        for field in fields:
            row.append(act_data[field])
        
        writer.writerow(row)

bioact_df_tmp = pd.read_csv('../data/T1/bioact.csv')

bioact_df = bioact_df_tmp.query('type == "IC50" and assay_type == "B" and relation == "="') # 条件に合致する行を残す
bioact_df.head(10)
activity_id assay_chembl_id assay_description assay_type molecule_chembl_id type units relation value target_chembl_id target_organism
0 32260 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL68920 IC50 uM = 0.041 CHEMBL203 Homo sapiens
1 32267 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL69960 IC50 uM = 0.170 CHEMBL203 Homo sapiens
2 32680 CHEMBL677833 In vitro inhibition of Epidermal growth factor... B CHEMBL137635 IC50 uM = 9.300 CHEMBL203 Homo sapiens
3 32770 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL306988 IC50 uM = 500.000 CHEMBL203 Homo sapiens
4 32772 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL66879 IC50 uM = 3000.000 CHEMBL203 Homo sapiens
5 32780 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL77085 IC50 uM = 96.000 CHEMBL203 Homo sapiens
6 33406 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL443268 IC50 uM = 5.310 CHEMBL203 Homo sapiens
7 34039 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL76979 IC50 uM = 264.000 CHEMBL203 Homo sapiens
8 34041 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL76589 IC50 uM = 0.125 CHEMBL203 Homo sapiens
9 34049 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL76904 IC50 uM = 35.000 CHEMBL203 Homo sapiens 
bioact_df.shape
# (7177, 11)

欠損値のあるエントリーを削除します。

bioact_df = bioact_df.dropna(axis=0, how = 'any')
bioact_df.shape
# (7176, 11)

重複を削除します。同じ化合物(molecule_chembl_id) について複数回の評価データが登録されていることが時々あります。今回は最初の結果だけを採用します。

bioact_df = bioact_df.drop_duplicates('molecule_chembl_id', keep = 'first')
bioact_df.shape
#  (5505, 11)

生理活性値がモル濃度の単位で評価されているものだけを残したいと思います。次のprint文は、どのような単位が含まれているかを確認し、行をいくつか削除したあとでどのようなエントリーが残すか調整するのに役立ちます。

print(bioact_df.units.unique())
bioact_df = bioact_df.drop(bioact_df.index[~bioact_df.units.str.contains('M')])
print(bioact_df.units.unique())
bioact_df.shape

# ['uM' 'nM' 'M' "10'1 ug/ml" 'ug ml-1' "10'-1microM" "10'1 uM"
#   "10'-1 ug/ml" "10'-2 ug/ml" "10'2 uM" '/uM' "10'-6g/ml" 'mM' 'umol/L'
#   'nmol/L']
#  ['uM' 'nM' 'M' "10'-1microM" "10'1 uM" "10'2 uM" '/uM' 'mM']
# (5435, 11)

訳注(04.2020時点)
上のセルの条件 ~bioact_df.units.str.contains('M')の~はnotを意味します。
pandasの条件ではand、or、notの代わりに&、|、~を使うことでエラーを回避できるそうです。参考
訳注ここまで

いくつか行を削除しましたが、あとでインデックスにしたがって反復処理（iteration）を行いたいので、インデックスを連続した値にリセットします。

bioact_df = bioact_df.reset_index(drop=True)  # drop=Trueで元の引数を削除
bioact_df.head()
activity_id assay_chembl_id assay_description assay_type molecule_chembl_id type units relation value target_chembl_id target_organism
0 32260 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL68920 IC50 uM = 0.041 CHEMBL203 Homo sapiens
1 32267 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL69960 IC50 uM = 0.170 CHEMBL203 Homo sapiens
2 32680 CHEMBL677833 In vitro inhibition of Epidermal growth factor... B CHEMBL137635 IC50 uM = 9.300 CHEMBL203 Homo sapiens
3 32770 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL306988 IC50 uM = 500.000 CHEMBL203 Homo sapiens
4 32772 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL66879 IC50 uM = 3000.000 CHEMBL203 Homo sapiens

IC50値をさらに比較可能とするために、全ての単位をnMに変換します。まず、ヘルパー関数を書きます。この関数を次工程でデータフレーム全体に適用します。

def convert_to_NM(unit, bioactivity):
#     c=0
# for i, unit in enumerate(bioact_df.units):
    if unit != "nM":        
        if unit == "pM":
            value = float(bioactivity)/1000
        elif unit == "10'-11M":
            value = float(bioactivity)/100
        elif unit == "10'-10M":
            value = float(bioactivity)/10
        elif unit == "10'-8M":
            value = float(bioactivity)*10
        elif unit == "10'-1microM" or unit == "10'-7M":
            value = float(bioactivity)*100
        elif unit == "uM" or unit == "/uM" or unit == "10'-6M":
            value = float(bioactivity)*1000
        elif unit == "10'1 uM":
            value = float(bioactivity)*10000
        elif unit == "10'2 uM":
            value = float(bioactivity)*100000
        elif unit == "mM":
            value = float(bioactivity)*1000000
        elif unit == "M":
            value = float(bioactivity)*1000000000
        else:
            print ('unit not recognized...', unit)
        return value
    else: return bioactivity

bioactivity_nM = []
for i, row in bioact_df.iterrows():
    bioact_nM = convert_to_NM(row['units'], row['value'])
    bioactivity_nM.append(bioact_nM)
bioact_df['value'] = bioactivity_nM
bioact_df['units'] = 'nM'
bioact_df.head()
activity_id assay_chembl_id assay_description assay_type molecule_chembl_id type units relation value target_chembl_id target_organism
0 32260 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL68920 IC50 nM = 41.0 CHEMBL203 Homo sapiens
1 32267 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL69960 IC50 nM = 170.0 CHEMBL203 Homo sapiens
2 32680 CHEMBL677833 In vitro inhibition of Epidermal growth factor... B CHEMBL137635 IC50 nM = 9300.0 CHEMBL203 Homo sapiens
3 32770 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL306988 IC50 nM = 500000.0 CHEMBL203 Homo sapiens
4 32772 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL66879 IC50 nM = 3000000.0 CHEMBL203 Homo sapiens

化合物データ

EGFRに対して評価を実施された全ての化合物（とそれぞれの評価方法）を含むデータフレームが手に入りました。それでは、それぞれのChEMBL IDの裏に隠れた化合物を取得したいと思います。

化合物のリストを取得する

ChEMBLから抽出し生理活性データを明確にした化合物を見てみましょう。まず、望みの生理活性データを有する化合物のChEMBL IDと構造を取り出します。

cmpd_id_list = list(bioact_df['molecule_chembl_id'])
compound_list = compounds.filter(molecule_chembl_id__in = cmpd_id_list)\
                              .only('molecule_chembl_id','molecule_structures')

次に、リストをPandasのDataFrameに変換し、重複を削除します。（くどいですが、padas from_records　関数の実行は時間がかかるかもしれません）。

訳注(04.2020時点)

ここでもクエリのデータ取得の段階での項目の絞り込み(only)がうまく機能しておらず、DataFrameへの変換が非常に遅いので、csvの出力ののちDataFrameに変換するという方法をとりました。出力したcsvは17MBあり1時間弱かかりました。

訳注ここまで

# 出力したいフィールド （csvのheaderにする）
cmpd_fields = ['molecule_chembl_id','molecule_structures']

# cmpd.csvというファイルを作成し出力する
with open('../data/T1/cmpd.csv', 'w') as f:
    writer = csv.writer(f)
    
    # ヘッダーの書き込み
    writer.writerow(cmpd_fields)
    
    # activityのデータを１行ずつ書き込み
    for cmpd in compound_list:
        row = []
        
        for field in cmpd_fields:
            row.append(cmpd[field])
        
        writer.writerow(row)

# compound_df = pd.DataFrame.from_records(compound_list)
compound_df = pd.read_csv('../data/T1/cmpd.csv')
compound_df = compound_df.drop_duplicates('molecule_chembl_id', keep = 'first')
print(compound_df.shape)
print(bioact_df.shape)
compound_df.head()
# (5435, 2)
# (5435, 11)
molecule_chembl_id molecule_structures
0 CHEMBL6246 {'canonical_smiles': 'O=c1oc2c(O)c(O)cc3c(=O)o...
1 CHEMBL10 {'canonical_smiles': 'C[S+]([O-])c1ccc(-c2nc(-...
2 CHEMBL6976 {'canonical_smiles': 'COc1cc2c(cc1OC)Nc1ncn(C)...
3 CHEMBL7002 {'canonical_smiles': 'CC1(COc2ccc(CC3SC(=O)NC3...
4 CHEMBL414013 {'canonical_smiles': 'COc1cc2c(cc1OC)Nc1ncnc(O...

この段階では、構造の情報について複数の異なる構造表現を保持しています。この中からカノニカルSMILESだけを残したいと思います。

訳注（04.2020時点）
オリジナルの方法とは異なりcsvからDataFrameにしたことで、カラム「molecule_structures」の要素は見た目はdict型ですがstr型となっています。
このままではcanonical_smilesをvalueとして取り出すことができません。
以下ではastモジュールのliteral_eval関数を適用してstr型からdict型に変換する処理を最初に実行しています。
また、forループで処理する際の、欠損値行のスキップがもともとのif条件ではうまくいかなかったため、float型の場合スキップするという条件に変更しています。
訳注ここまで

# mapで型変換をデータフレーム全体に適用
# 欠損値でエラーが出るためmapの引数 na_actionを利用

import ast
compound_df['molecule_structures'] = compound_df['molecule_structures'].map(ast.literal_eval, na_action='ignore') 

for i, cmpd in compound_df.iterrows():
    # オリジナルのif条件
    # if compound_df.loc[i]['molecule_structures'] != None:
    
    if type(compound_df.loc[i]['molecule_structures']) != float: # Nanがfloat型であることを利用したif条件に書き換えた
        compound_df.loc[i]['molecule_structures'] = cmpd['molecule_structures']['canonical_smiles']

print (compound_df.shape)
# (5435, 2)

出力データの準備

知りたい値をChEMBL IDに基づいて一つのデータフレームに統合します。 * ChEMBL-IDs * SMILES * 単位 * IC50

output_df = pd.merge(bioact_df[['molecule_chembl_id','units','value']], compound_df, on='molecule_chembl_id')
print(output_df.shape)
output_df.head()
# (5435, 4)
molecule_chembl_id units value molecule_structures
0 CHEMBL68920 nM 41.0 Cc1cc(C)c(/C=C2\C(=O)Nc3ncnc(Nc4ccc(F)c(Cl)c4)...
1 CHEMBL69960 nM 170.0 Cc1cc(C(=O)N2CCOCC2)[nH]c1/C=C1\C(=O)Nc2ncnc(N...
2 CHEMBL137635 nM 9300.0 CN(c1ccccc1)c1ncnc2ccc(N/N=N/Cc3ccccn3)cc12
3 CHEMBL306988 nM 500000.0 CC(=C(C#N)C#N)c1ccc(NC(=O)CCC(=O)O)cc1
4 CHEMBL66879 nM 3000000.0 O=C(O)/C=C/c1ccc(O)cc1

カラム（列）名を、それぞれIC50とSMILESに変更します。

output_df = output_df.rename(columns= {'molecule_structures':'smiles', 'value':'IC50'})
output_df.shape
# (5435, 4)

SMILES表記が手に入らない化合物については、以降のトークトリアルで使用することができません。ですので、SMILES列に値のない化合物を削除します。

output_df = output_df[~output_df['smiles'].isnull()]
print(output_df.shape)
output_df.head()
# (5428, 4)
molecule_chembl_id units IC50 smiles
0 CHEMBL68920 nM 41.0 Cc1cc(C)c(/C=C2\C(=O)Nc3ncnc(Nc4ccc(F)c(Cl)c4)...
1 CHEMBL69960 nM 170.0 Cc1cc(C(=O)N2CCOCC2)[nH]c1/C=C1\C(=O)Nc2ncnc(N...
2 CHEMBL137635 nM 9300.0 CN(c1ccccc1)c1ncnc2ccc(N/N=N/Cc3ccccn3)cc12
3 CHEMBL306988 nM 500000.0 CC(=C(C#N)C#N)c1ccc(NC(=O)CCC(=O)O)cc1
4 CHEMBL66879 nM 3000000.0 O=C(O)/C=C/c1ccc(O)cc1

次のセルを実行すると、小さなIC50値が読みづらいことがわかると思います。そこで、IC50値をpIC50値に変換したいと思います。

output_df = output_df.reset_index(drop=True)
ic50 = output_df.IC50.astype(float) 
print(len(ic50))
print(ic50.head(10))

5428
0         41.0
1        170.0
2       9300.0
3     500000.0
4    3000000.0
5      96000.0
6       5310.0
7     264000.0
8        125.0
9      35000.0
Name: IC50, dtype: float64

# IC50をpIC50に変換し、pIC50の列を付け足す
pIC50 = pd.Series() 
i = 0
while i < len(output_df.IC50):
    value = 9 - math.log10(ic50[i]) # pIC50=-log10(IC50 mol/l) --> for nM: -log10(IC50*10**-9)= 9-log10(IC50)
    if value < 0:
        print("Negative pIC50 value at index"+str(i))
    pIC50.at[i] = value
    i += 1
    
output_df['pIC50'] = pIC50
output_df.head()
molecule_chembl_id units IC50 smiles pIC50
0 CHEMBL68920 nM 41.0 Cc1cc(C)c(/C=C2\C(=O)Nc3ncnc(Nc4ccc(F)c(Cl)c4)... 7.387216
1 CHEMBL69960 nM 170.0 Cc1cc(C(=O)N2CCOCC2)[nH]c1/C=C1\C(=O)Nc2ncnc(N... 6.769551
2 CHEMBL137635 nM 9300.0 CN(c1ccccc1)c1ncnc2ccc(N/N=N/Cc3ccccn3)cc12 5.031517
3 CHEMBL306988 nM 500000.0 CC(=C(C#N)C#N)c1ccc(NC(=O)CCC(=O)O)cc1 3.301030
4 CHEMBL66879 nM 3000000.0 O=C(O)/C=C/c1ccc(O)cc1 2.522879

EGFRを標的として集めた生理活性データ

集めたデータセットを眺めてみましょう。

化合物を描画する

次の工程ではデータフレームに化合物を追加し、もっとも大きいpIC50値をもつ化合物の構造を眺めます。

PandasTools.AddMoleculeColumnToFrame(output_df, smilesCol='smiles')

pIC50値で化合物をソートします。

output_df.sort_values(by="pIC50", ascending=False, inplace=True)
output_df.reset_index(drop=True, inplace=True)

もっとも活性な分子を、もっとも大きなpIC50値の分子として描画します。

output_df.drop("smiles", axis=1).head()
molecule_chembl_id units IC50 pIC50 ROMol
0 CHEMBL63786 nM 0.003 11.522879 Mol
1 CHEMBL53711 nM 0.006 11.221849 Mol
2 CHEMBL35820 nM 0.006 11.221849 Mol
3 CHEMBL53753 nM 0.008 11.096910 Mol
4 CHEMBL66031 nM 0.008 11.096910 Mol

訳注(04.2020時点)

ROMolに構造が描画されず、（SMILESでもない）文字列になっている場合はRDKitとPandasのバージョンの互換性がよくない可能性があります。Pandasを version 0.25.0 以前のものにすることで解決する可能性があります。 RDKit イシュー　を参照してください。

訳注ここまで

出力ファイルの書き出し

以降のトークトリアルでデータを使うため、csvファイルとしてデータを保存します。データを保存するときには化合物の列（分子の構造図のみをもつ列）を削除することをお勧めします。

output_df.drop("ROMol", axis=1).to_csv("../data/T1/EGFR_compounds.csv")

ディスカッション

このチュートリアルではChEMBLのデータベースから、私たちの目的とするターゲット分子に対して取得可能な全ての生理活性データを集めました。データセットをIC50あるいはpIC50の測定値をもつものだけにフィルタリングしました。

ChEMBLのデータはさまざまなデータソースに由来するものであることに注意してください。化合物のデータは世界中の様々なラボの様々な人々によって生み出されたものです。したがって、このデータセットを使って予測を行うときには注意する必要があります。結果を解釈し、自分たちの予測にどれだけの信頼を置けるか決める際には、データのソースとデータの生み出されたアッセイの一貫性について、常に慎重に考えることが重要です。

次のチュートリアルでは、手に入れたデータをリピンスキーのルールオブファイブでフィルタリングし、そして好ましくない部分構造でフィルタリングします。もう一つな他の重要なステップとして、データを綺麗に整型し、重複を取り除くことが挙げられます。（まだ、）私たちのトークトリアルでは取り上げられていないので、このタスクに使えるツールとして、標準化ライブラリ (github Francis Atkinson) あるいは MolVS というものがあることに言及しておきたいと思います。

クイズ

  • このトークトリアルではChEMBLから化合物と生理活性データをダウンロードしました。ChEMBLデータベースは他にどのように利用することができるしょうか？
  • IC50とEC50の違いは何でしょうか？
  • ChEMBLから取り出したデータはどのような目的で使いことができるでしょうか？

訳以上

追記

私はCADD、プログラミングの専門性がないので翻訳、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。 また環境の違いのためか、GitHubのコードそのままでは実行できない箇所が散見されました。行き当たりばったりでごり押ししているので誤り等多々あると思います。 オススメの解決策等あればご教示いただければ幸いです。

Jupyter notebookからmarkdwonで出力ごとダウンロードできるの知らなかった！便利！

はじめに ~TeachOpenCADD 環境構築~

TeachOpenCADD

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

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

github.com

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

以下、訳。

TeachOpenCADD

TeachOpenCADD: オープンソースのパッケージとデータを使ったコンピューター支援（CADD）ドラッグデザインの教育用プラットフォーム

Dominique Sydow, Andrea Morger, Maximilian Driller and Andrea Volkamer

In Silico Toxicology, Institute for Physiology, Universitätsmedizin Berlin, Virchowweg 6, 10117 Berlin

DOI DOI Binder

質問がある場合は teachopencadd@charite.de に連絡してください。もしくは、このページにない講義用チュートリアルについての情報、チュートリアルのもととなるアイデアがあればぜひ連絡してください。

このページは　CC BY 4.0 ライセンスののもとで提供されています。ライセンスの詳細は　http://creativecommons.org/licenses/by/4.0/　をご覧ください。

目次

  • この教育用プラットフォームの目的
  • トピック
  • 使用方法
  • 連絡先
  • ライセンス
  • 引用

この教育用プラットフォームの目的

コンピューター支援ドラッグデザイン (CADD)にとって、ケモインフォマティクスと構造バイオインフォマティクスに関するオープンソースのプログラミングパッケージは、パワフルなツールで、モジュール式かつ再現可能、簡単に共有可能なパイプラインを提供します。ところが、ツールの説明、ドキュメンテーションは利用可能な一方で、CADDの活用方法について焦点をあて、背景にあるコンセプトを教える上で、自由に利用可能な具体例はほんの少ししかありません。例を挙げれば TDT initiative と言ったものがあります。 [1]

ここでは、学生による学生のために開発されたCADDの教育用プラットフォームを提供します。使用するオープンソースパッケージはPythonツール、RDKit [2] 、PyPDB [3] 、そして PyMol [4] です。それぞれのトピックについて、インタラクティブなJupyter Notebook [5] を開発しました。Githubで自由に利用可能で、トピックに関する理論的な背景の詳細な説明と、詳しいコメント付きのPythonコードが書いてあります。EGFR キナーゼを例に、公共の化合物データベース ChEMBL [6] からデータを取得する方法、ドラッグライクネスの観点から化合物をフィルタリングする方法、好ましく無い部分構造をもつ化合物を探し出し取り除く方法を議論します。さらに、化合物の類似性を評価する方法を導入し、化合物のクラスター化を行い、多様性のある化合物セットをつくります。また、新規な活性化合物を予測するためのモデルを構築するため、機械学習のアプローチを導入します。最後に、公共のタンパク質データベース PDB [7] から構造を取得し、それらを使ってリガンドベースの3D ファーマコフォアの生成、オフターゲットを予測するために結合サイト比較を実施します。

このプラットフォームを作成したことの目的は、オープンソースのケモインフォマティクスツールに興味を持っている利用者に、いかに簡単に利用できて、どんな利点があるかということを紹介することです。取り上げるトピックは継続的に増やしていきますし、コミュニティからのコントリビューションもお待ちしています。教育目的の利用だけでなく、この Jupiter notebookは各プロジェクトに合わせて修正し拡張するための出発点としても使えます。

文献:

[1] S. Riniker et al., F1000Research, 2017, 6, 1136 [2] G. Landrum, RDKit [3] W. Gilpin, Bioinformatics, 2016, 32, 156-60 [4] The PyMOL Molecular Graphics System, Version 1.8, Schrödinger, LLC [5] T. Kluyver et al., IOS Press, 2016, 87-90. [6] A. Gaulton et al., Nucleic Acid Res., 2017, 40, D1100-7 [7] H. Berman et al., Nucleic Acid Res., 2000, 28, 235-42

トピック

f:id:magattaca:20200405221750j:plain
TeachOpenCADD文献（D. Sydow et al., J Chem, 2019, 11, 29）Figure 1より

TeachOpenCADDのトピックはいわゆる トークトリアル（talkatorial） の形式（理論とコードのミックスで、プレゼンテーション形式での発表にも利用可能な形式）で取り上げられていて、現在、以下の内容で構成されています：

ケモインフォマティクス :

  • Topic 1. 化合物データの取得：ChEMBL
  • Topic 2. 化合物フィルタリング：ADMEとリードライクネスのクライテリア
  • Topic 3. 化合物フィルタリング：好ましく無い部分構造
  • Topic 4. リガンドベーススクリーニング: 化合物類似性
  • Topic 5. 化合物クラスタリング
  • Topic 6. 最大共通部分構造（Maximum common substructures）
  • Topic 7. リガンドベーススクリーニング：機械学習

構造バイオインフォマティクス：

  • Topic 8. タンパク質データの取得：Protein Data Bank (PDB)
  • Topic 9. リガンドベースファーマコフォア
  • Topic 10. 結合サイトの類似性
  • Topic 11. オンライン API/サーバー を使用した構造ベースのCADD
    • 11a. 潜在的なキナーゼ阻害剤を見つけるためのKLIFSとPubChemのクエリ化
    • 11b. 標的に対して候補のドッキングを行う
    • 11c. 結果の可視化と既知のデータとの比較

教育用の素材は次のフォーマットで提供されています：

  • コーディングベースのJupyter Notebook（topics 1-11）はGithubのこのページにあり、いわゆるトークトリアル（talktorial = talk + tutorial）、つまりプレゼンテーションでも使えるチュートリアルの形式です。
  • GUIベースのKNIME ワークフロー (topics 1-8) は KNIME Hub で手に入ります。

使用方法

レポジトリと必要な依存関係にあるものを全てホストするためにBinderを使うことができます。あるいはローカルでトークトリアルを使うこともできます（レポジトリをダウンロートして、依存関係にあるものをインストールしてください）

Binderの使い方 (対象トークトリアル T1-T8)

以下のリンクに従うことでレポジトリと必要な依存関係にあるものを全てホストするためにBinderを使うことができます:

Binder

セットアップには２、３分かかります。

ローカル環境へのインストール (対象トークトリアル　全て)

Linux

  • 1（トークトリアルを含む）TeachOpenCADDレポジトリのコピーをローカルに落としてください

    1. ... zipアーカイブとしてダウンロードしてunzipしてください： f:id:magattaca:20200405221951p:plain

    2. ... もしくは git パッケージを使って自分のコンピュータにクローンしてください：

git clone https://github.com/volkamerlab/TeachOpenCADD.git

  • 2 Anacondaソフトウェアを使ってパッケージのバージョン管理をクリーンにしてください。

    Anaconda2、あるいはAnaconda3をインストールします。理論上はインストールするAnacondaのバージョンの問題です。ですが、Anaconda2でしかテストしていないので、Anaconda3で同じような動作となるかは保証できません。

    https://docs.anaconda.com/anaconda/install/

  • 3 パッケージ管理システムcondaを使ってトーカトリアルの環境を作ります（ teachopencadd という名前をつけています）

    必要なパッケージをすべて含んでいる環境ファイル (yml file) を提供しています。

conda env create -f environment.yml

Note: この環境を手動でつくることもできます。 "別の方法：conda環境を手動で作る" を参照してください。

  • 4 作成したconda環境をアクティベートしてください。
conda activate teachopencadd

これでconda環境の中で作業することができます。

  • 5 作成したconda環境をJupyter notebook にリンクしてください
python -m ipykernel install --user --name teachopencadd

# FYI: このリンクをアンインストールするには以下のコマンドを使用してください:
#jupyter kernelspec uninstall teachopencadd
  • 6 Jupyter notebookを起動してください。
jupyter notebook

  • 7 メニューからJupyterのカーネルを作成したconda環境に変えてください。

    Kernel > Change kernel > Python [conda env:teachopencadd]

    f:id:magattaca:20200405222419p:plain

  • 8 これで最初のトーカトリアル開始です。Enjoy!

別の方法：conda環境を手動で作る

Note: これは上のステップ３で書かれているymlファイルを使用してconda環境をつくる方法の代わりとなる方法です。

# `teachopencadd` という名前の環境をつくりアクティベートします
conda create -n teachopencadd python=3.6
conda activate teachopencadd

# conda経由でパッケージをインストールします
conda install jupyter  # ipykernelもインストールします
conda install -c conda-forge rdkit  # umpy と pandas もインストールします
conda install -c samoturk pymol  # Linux: freeglut と glewもインストールします
conda install -c conda-forge pmw  # Pymolのターミナルウィンドウを立ち上げるのに必要です。
conda install -c conda-forge scikit-learn  # scipyをインストールします
conda install -c conda-forge seaborn  # matplotlibをインストールします
conda install -c conda-forge chembl_webresource_client
conda install -c conda-forge biopandas
conda install -c conda-forge pypdb
#conda install jupyter ipykernel pandas scikit-learn seaborn  # おそらくすでにインストールされていると思います。

# pip経由でパッケージをインストールします (おそらくデフォルトインストールされています）

pip install chembl_webresource_client biopandas pypdb

Windows

大枠では、LinuxとWindowsで同じように問題なくインストールできます。

唯一の違いはPyMolのインストール作業です。

cd C:/Anaconda2
set path=%path%;C:\Anaconda2
pip install pymol‑2.3.0a0‑cp36‑cp36m‑win_amd64.whl
pip install pymol_launcher‑2.1‑cp36‑cp36m‑win_amd64.whl
# 上の文言はダウンロードした pymol_XXXX.whl ファイルに合わせてください

MacOS

Linux環境と同様にインストールできますが、pymol をオープンソース samoturk condaチャンネルからインストールすることはできませんでした。代わりにschrodinger チャンネルを使うことができます。残念ながらPyMOLを実行するにはSchrödingerライセンスが必要です（ライセンスは 教育目的 の使用には無料で提供されています）。

インストールを手動で行うには、次のコマンドを：

conda install -c samoturk pymol # Installs also freeglut and glew

次に置き換えてください

conda install -c schrodinger pymol  # Installs also freeglut and glew

出くわした問題

RDKit
  • from rdkit.Chem.Draw import IPythonConsole がエラーメッセージ ImportError: No module named 'PIL'を投げる: pip install pillow のインストールを試してください(https://github.com/rdkit/rdkit/issues/1179 を確認してみてください)
  • Draw.MolsToGridImage が Serial errorをだす: データのタイプを Draw.MolsToGridImage(list(df.ROMol), useSVG=True) のようにリストしてください。
pip

バージョン　10.0.X 以降のパッケージマネージャーは以前のものと動作が異なっています。上で使われているコマンドのpip install chembl_webresource_clientは importエラー cannot import name mainを返します:

この問題はpython経由でpipを呼び出すことで解決することができます：

python3 -m pip install chembl_webresource_client
python3 -m pip install pypdb

chembl_webresource_client をインストールする際に、 urllib3 の互換性がないという警告が出るかもしれません。アップデートすることで解決することができます：

python3 -m pip install urllib3 --upgrade

解決方法は次からとりました。:

https://stackoverflow.com/questions/28210269/importerror-cannot-import-name-main-when-running-pip-version-command-in-window

連絡先

質問や提案があれば連絡してください！

皆様からのご意見をお待ちしています！

ライセンス

この仕事は Attribution 4.0 International (CC BY 4.0) のもとでライセンスされています。このライセンスのコピーをみるには次のリンクを参照してください http://creativecommons.org/licenses/by/4.0/ 。

引用

TeachOpenCADDプラットフォームの著者らは次の出資者からの公的出資を受け取っています:

  • Bundesministerium für Bildung und Forschung (Grant Number 031A262C)
  • Deutsche Forschungsgemeinschaft (Grant number VO 2353/1-1)
  • HaVo-Stiftung, Ludwigshafen, Germany
  • Stiftung Charité (Einstein BIH Visiting Fellow project)
  • "SUPPORT für die Lehre" program (Förderung innovativer Lehrvorhaben) of the Freie Universität Berlin
  • Open Access Publication Fund of Charité – Universitätsmedizin Berlin

サイエンティフィックな文献で TeachOpenCADD の教材を利用した時は、我々の文献を引用してください：

TeachOpenCADDプラットフォームの有用さを計り、今後の出資をうけるために役立ちます！

@article{TeachOpenCADD2019, author = {Sydow, Dominique and Morger, Andrea and Driller, Maximilian and Volkamer, Andrea}, doi = {10.1186/s13321-019-0351-x}, journal = {J. Cheminform.}, number = {1}, pages = {29}, title = {{TeachOpenCADD: a teaching platform for computer-aided drug design using open source packages and data}}, url = {https://doi.org/10.1186/s13321-019-0351-x}, volume = {11}, year = {2019} }

@article{TeachOpenCADDKNIME2019, author = {Sydow, Dominique and Wichmann, Michele and Rodr{\'{i}}guez-Guerra, Jaime and Goldmann, Daria and Landrum, Gregory and Volkamer, Andrea}, doi = {10.1021/ACS.JCIM.9B00662}, journal = {J. Chem. Inf. Model.}, publisher = {American Chemical Society}, title = {{TeachOpenCADD-KNIME: A Teaching Platform for Computer-Aided Drug Design Using KNIME Workflows}}, url = {https://pubs.acs.org/doi/full/10.1021/acs.jcim.9b00662}, year = {2019} }

訳以上

追記

以上でREADMEの訳となります。環境構築に関して私の場合、MacBook Pro(2017) macOS High Sierraを使用しており、conda環境での構築を試しています。
他のOSでのインストールやyml、Binder（って何かわかりませんでした・・・）を試していないので、不具合の有無はわかりません。
ソフトウェアのバージョンアップはそれぞれまちまちなので、場合によっては互換性がなくなってしまっているものもあるかもしれません。
ぜひぜひ情報をシェアしてください。（…このページについて本家の連絡先には問い合わせないでください）

素人の試訳なので誤訳、間違い等あると思います。ご指摘いただければ幸いです。

・・・いつかGitHubにちゃんとjupyter notebookをあげたい。

オープンソースで医薬品探索のコンセプトを勉強しよう！〜TeachOpenCADD〜

TeachOpenCADD

TeachOpenCADDをご存知でしょうか？ざっくり言うとCADDの自習用教材です。*1

f:id:magattaca:20200405201048p:plain:w300
TeachOpenCADD

CADD（Coputer-Aided Drug Design、コンピュータ支援医薬品設計）というのは、 その名の通り「医薬品の探索、ドラッグデザインにコンピュータを利用する」ということです。（≒ in silico創薬？）

医薬品探索における計算機利用、シミュレーションの重要性とともにその門戸はどんどんと開かれており、
1. 計算機性能の向上
2. 利用可能な（オープンソース）ソフトウェアが増えた
3. 公開データベース
などがその理由としてよく挙げられています。

それぞれ説明書もあって道具は揃った！よし、あとは実際につかうだけ！・・・とはいかないですよね。
とりわけ私のような初学者にとっては「車あるけどカーナビもスマホも無い。どっち行けばいいの？？」状態です。

ツールと実践の間のハードルをクリアするためには、
1. そもそもなぜそのツールが必要とされたのか、背景にあるコンセプトを学ぶこと
2. 実践的な具体例でツールを利用してみること
が必要です。

つまり、あれです。「先達・・・あらまほし・・・」ってやつですね。

で、そんな迷える仁和寺の法師たちのためのプロジェクトがこちらTeachOpenCADD platformです。

TeachOpenCADD?

TeachOpenCADDの目的は「CADDに馴染みのない学生、研究者のために、そのコンセプトと実践を学習するための教材を提供する」ことで、
チュートリアル形式で一歩ずつ学習を進められるように設計された、学生による学生のためのプラットフォームとなっています。

開発の主体となっているのはシャリテ（The Charité、ベルリン大学病院）のVolkamer研究室で、 構造ベースの手法による毒性予測やオフターゲット予測を主な研究テーマとされているようです。

TeachOpenCADDの特徴として、

  • 無料で利用可能なソースを使用
  • Interactiveに学習できるJupyter notebook形式 (GitHubに公開)

となっており、「EGFRキナーゼ阻害剤の探索」を具体例にしてリガンドベース、構造データベースのCADDの手法を順に学ぶことができるようになっています。

f:id:magattaca:20200405202008p:plain
TeachOpenCADDで学べること（J Cheminform(2019)11:29 Fig1.より）

オープン（Open）ソースでCADDを教える(Teach)教材、TeachOpenCADDということのようですね。

各チュートリアルは自習だけでなく講義で説明しながら使えるようにもなっており、 トークトリアル（talk + tutorial = talktorial)とも呼ばれています。

f:id:magattaca:20200405202310p:plain
実際のチュートリアル（J Cheminform (2019)11:29 Fig2）

以上は私の雰囲気翻訳なので正確な情報はオリジナルを参照してください。

  1. TeachOpenCADD Jupyter Notebooks: Talktorials 1-10.
  2. TeachOpenCADD KNIME workflows
  3. GitHub

EGFRキナーゼ阻害剤

さて、TeachOpenCADDではチュートリアル 10回（GitHubの最新版は11回まで）を通してCADDのコンセプトを学びますが、 一貫して取り上げられている題材「EGFRキナーゼ阻害剤（EGFR-TKI）」というのはどういうものなのでしょうか？

まずは敵を知る必要がある！！*2

Wikipediaによると、EGFR（上皮成長因子受容体、Her1）は細胞膜上の受容体で、膜貫通型糖タンパク質です。細胞内にチロシンキナーゼ活性を有するドメインがあり、細胞外からの刺激を細胞内でリン酸化シグナルとして伝達します。このシグナルは細胞の分化・増殖に関わっており、したがって機能変異により発がん、がんの増殖・浸潤・転移に関与するようになります。*3
EGFR遺伝子の変異は肺がんで多く見られるドライバー変異で、日本人の場合非小細胞肺がん患者で多く変異が見られているそうです。*4

変異による過剰な増殖シグナルをとめることで抗がん作用が期待できるということで、細胞外（抗EGFR抗体）および細胞内（EGFRキナーゼ阻害剤）の２つのアプローチで医薬品の開発、上市されています。

EGFRキナーゼ阻害剤といえばやはり一番最初に思い当たるのはゲフィチニブではないしょうか？薬学においてゲフィチニブは様々な点で取り上げられる医薬品です。

  • 初期の分子標的薬
  • 優先審査により世界に先駆けて日本で承認
  • 副作用（イレッサをめぐる問題）
  • 遺伝子変異と薬効の関係

など様々な文脈で名前を聞いたように思います。*5

すでに第3世代までが承認されているEGFRキナーゼ阻害剤ですが、耐性発現、既存の薬剤で効果のでない変異があるなど依然としてアンメットニーズがあるそうです。*6 また、EGFR陽性の肺がんでは免疫チェックポイント阻害薬の効果が出にくいという報告もあるそうで、 EGFR阻害剤最初の承認から20年近くたった今でもまだまだ研究する余地がありそうです。*7

で、私は何がしたいのか？

Q: ぐだぐだと書いているが結局何がしたいのか？
A: 一人で勉強するのが寂しい

ということでのんびり翻訳して晒していこうと思います。 ライセンスCC BY 4.0だから怒られないはず・・・

ではでは！

~~ (追記 2020/09/28) ~~~

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

github.com

*1:文献. TeachOpenCADD Jupyter Notebooks: Talktorials 1-10.

*2:私は医療従事者でも専門家でもないので記述の正確性に問題があることをご承知おきください。

*3:上皮成長因子受容体_Wikipedia

*4:日本内科学会雑誌2014(103)1355

*5:このあたり付け焼き刃の知識で踏み込めないので興味を持たれた方は調べてみてください。

*6:【肺がん】EGFR-TKI 一次治療に相次ぎ新薬…それでも残るアンメットニーズ

*7:名古屋大学プレスリリース

pymolをimportしようとする話の続きの話

PyMOLのjupyter notebookからのimportに失敗したことについて記事を書いたところ、早速twitterで解決方法を教えていただきました。 @kira_kira_worldさん、@Ag_smith先生ありがとうございました！

解決策に入る前に・・・pythonのバージョンについて

先の記事を書いた時点では気づいていなかったのですがどうもpythonのバージョンに原因があったようです。

私の環境はmacOS High Sierraで

  • Homebrew > pyenv > Anaconda（複数の仮想環境）

みたいな感じになっています。（たぶん）

で、全体（pyenv global）はPython 3.70となっています。 この状態でhomebrewからPyMOLをインストールしたためか、importしようとしたPyMOLは以下のディレクトリ 「/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages」 を呼びに行っていました（先の記事でpathを通した所）。
如何にもpython3.7に紐づいていそうな場所です。

また、別の話として、色々とソフトウェアをインストールするためにAnacondaで複数の仮想環境を作成し、 jupyter notebookのカーネルに追加、切り替えつつ使おうとしています。どうやらここに問題がありそうでした。

先の記事で使用した仮想環境（カーネル）のpythonのバージョンを確認したところPython 3.6.10となっていました。 そこで別のPython 3.70の仮想環境に切り替えたところ、エラーは出ずにimport pymolすることができました。（pathを通す必要はありました）

おそらくPythonのバージョンをきちんと合わせておけば問題なく使用できるのではないでしょうか？ また確認不足で誤情報を垂れ流してしまった。。。すみません。

解決策1）シンボリックリンクの作成

ではでは、教えていただいた方法をご紹介したいと思います！
まずは@Ag_smith先生に教えていただいた方法です。

先の記事では怪しげなファイル「_cmd.cpython-37m-darwin.so」をコピペして「_cmd.so」に名前を変えるという荒技を行いました。 これだと思わぬ不具合を起こす可能性があるということで、こういう場合はシンボリックリンクを作成すれば良いそうです。

シンボリックリンクはUNIX系のOS（MacやLinux）で、Windowsのショートカットに相当するものだそうです。*1*2

つまり「_cmd.cpython-37m-darwin.so」のシンボリックリンクとして「_cmd.so」を作成しておけば、 import pymolはシンボリックリンクを中継として本体のファイルに辿りつくことができる、ということのようです。

作成方法は簡単で、ターミナルで以下を実行すればOKです。

# シンボリックリンクを作成したいファイルが入っているディレクトリに移動
cd /usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages/pymol

# 「ln -s 本体のファイル　リンクの名前」でシンボリックリンクを作成
ln -s _cmd.cpython-37m-darwin.so _cmd.so

できました！

f:id:magattaca:20200321225101p:plain

だそうです。「種類：エイリアス」となっていますがエイリアスとシンボリックリンクは微妙に機能が違うらしい・・・。 このあとjupyter notebookで（エラーの出ていたpython 3.6の仮想環境でも）import pymolできることを確認しました。

スマート！

解決策2）インストール方法を変える（Anaconda CloudのSchrödinger社のチャンネルを利用）

さて、もう一つの解決方法は@kira_kira_worldさんに教えていただいたもので、 condaでSchrödinger社のチャンネルからPyMOLをインストールするというものです。

１行でOK!

conda install -c schrodinger pymol

ですが注意点があります。こちらはSchrödinger社が配布しているものだということです。

Wikipediaによると、 もともとPyMOLはWarren L. DeLano氏により開発されたオープンソースのソフトウェアで、 現在はSchrödinger社によって商品化されているとのことです。 ですので、上記の方法で入手可能なPyMOLはライセンスに注意が必要そうです。

私はライセンス等に詳しくなく判断ができないので、以下手に入る情報を共有させていただきます。（←コンプライアンスを意識した発言）

まずはSchrödinger社のPyMOLのページ。
2020年03月20日現在、こちらではPyMOL 2.3のダウンロードが可能です（Python 3.7を含むそう）。 Windows、macOS、Linuxそれぞれのインストーラ、および上記のAnaconda Channelを利用する方法が提供されています。 また、ライセンスの購入方法についてもリンクがあります。

さらにスクロールしていくとOpen-Source Philosophyの文字が・・・
Schrödinger社が商用として維持、開発しているものの、元々のDelano氏の意向を引き続き大枠をオープンソースプロジェクトとしているようです。 したがってほとんどのソースコードはGitHub上で公開されています。ページのリンクはこちら。

Homebrew経由でインストールを行った際にも、ソースコードは上記のGithubから取得しているようです。 @Ag_smith先生の記事に記載されていて、 GitHubのHomebrewのFomulraでも取得先URLを確認できます。

では、Schrödinger社コンパイル済みのものをインストールした場合どうなるのでしょうか？試してみます。

テストのための仮想環境を作成し、condaで取得します。ちなみにAnaconda CloudのPyMolページはこちら。

# pythonのバージョンを指定して仮想環境pymoltestを作成
conda create -n pymoltest python=3.7

# 仮想環境をactivate
conda activate pymoltest

# PyMOLをインストール
conda install -c schrodinger pymol

# PyMOLを起動
pymol

無事PyMOLのウインドウが立ち上がりました！

f:id:magattaca:20200321225337p:plain

ただしHomebrew経由のインストールとは明確な違いが！

同時にActivationのウィンドウが立ち上がります。
「オフィシャル版なのでライセンスを持っているならactivationしてくださいね」ということのようです。 ライセンスがない場合、評価版として30日間使用することができるようです（skip activation）。

評価期間がすぎると「ウィンドウにNo License Fileの表示がされ使えなくなる機能がある」、とのことですが もうすでにNo License Fileの文字が見えますし、30日すぎるとどうなるのかはまだわかりません。

ちなみに学生であればライセンスは無料で取得することができるそうです。が、私はとうに学生ではなくなってしまったのでとりあえずskip activationで先に進めたいと思います。

さて、こちらのPyMOLがjupyter notebookでimportできるか検証したいと思います。 まず準備としてcondaの仮想環境をjupyter notebookに認識させます。

# 以下ターミナルで実行
# jupyter notebookのカーネルに仮想環境を追加
python -m ipykernel install --user --name pymoltest

# 一旦、仮想環境を抜けます
conda deactivate

また、インストールしたPyMOLの場所を確認するため、PyMOLのコマンド枠で以下を打ち込みます。

import sys
print(sys.path)

今回、それっぽいpathは「'/Users/XXXXX/.pyenv/versions/anaconda3-5.3.1/envs/pymoltest/lib/python3.7/site-packages'」でした（XXXXXは実際にはユーザーアカウント名ですが恥ずかしいので隠しています）。 condaの仮想環境pymoltestの中に入っていることがわかりやすいですね。

準備できたのでjupyter notebookを起動、新しいnotebookを作成してChange kernelでkernelをpymoltestに切り替えます。

f:id:magattaca:20200321225534p:plain

前回同様、pathを通してpymolをimportしてみます。

import sys
sys.path.append('/Users/XXXXX/.pyenv/versions/anaconda3-5.3.1/envs/pymoltest/lib/python3.7/site-packages') # XXXXXは実際にはユーザーアカウント名

import pymol

# importしたpymolのpathを確認
print(pymol.__path__)

無事importすることができました!!　PATHも指定したものとなっていました！

ちなみにpython 3.6の仮想環境のkernelで上記を実行したところ、前回記事同様に「モジュールpymol._cmdが見つからない」というエラーがでてimportできませんでした。 やはりバージョンの確認が大事そうですね。

何はともあれconda経由でもPyMolのインストールとCUIの利用が可能であることが確認できました！

目的は達成できたので、とりあえずオープンソース版での利用することとして仮想環境ごと削除しておきます。

# 以下ターミナルで実行
# jupyter notebookのカーネルからアンインストールしておく
jupyter kernelspec uninstall pymoltest

# condaの仮想環境を削除
conda remove -n pymoltest --all

ターミナルを一度落としてから再度立ち上げpymolを実行することでオープンソース版のPyMOLが立ち上がることを確認できました。

_cmdはどこからきたのか?

さて、jupyter notebookからPyMOLをimportする際の不具合は

  • Pythonのversion違いである可能性がある

こと、解決策として

があることまでわかりました。

で、今回importエラーとなった「pymol._cmd」とは何者か？もう少し調べてみたいと思います。

Github Open-Source PyMOL上でオープンソースのソースコードを確認してみましたが、 ディレクトリmodules/pymolの中にcmd.pyはあるものの_cmdに相当するものはありませんでした。

ビルド（？）する際に作成されたのかな？ということでsetup.pyの中身を確認したところ以下のコードに行き当たりました。

from distutils.core import setup, Extension
# ~ 中略 ~
ext_modules += [
    Extension("pymol._cmd",
              get_sources(pymol_src_dirs),
              include_dirs = inc_dirs,
              libraries = libs,
              library_dirs = lib_dirs,
              define_macros = def_macros,
              extra_link_args = ext_link_args,
              extra_compile_args = ext_comp_args,
              extra_objects = ext_objects,
    ),
# ~以下略~

distutils はPythonのモジュール配布のためのユーティリティで *3、Extensionは「C/C++拡張モジュールをセットアップスクリプトで表すために使われ」るそうです*4。

Extensionの説明と引数(argument)の内容からソースコードからビルドに伴って作成されたC/C++関連のファイルということがなんとなくぼんやりとわかりました。 これを踏まえて「_cmd.cpython-37m-darwin.so」というファイル名を再度眺めてみます。

前回邪魔なのではと削除してしまった「cpython-37m-darwin」の部分ですが、こちらはCythonでコンパイルする際に追加される接尾辞(suffix)のようです*5。 CythonはPythonの「処理速度が遅い」という問題を解決するためにC/C++に変換するもの（？？？？）で*6、 また、darwinはmacOS（Mac OS X）のベースとなっているオープンソースのオペレーティングシステムDarwinのことのようです*7。

なんとなくですがMacでpython 3.7でC/C++関連のコンパイルを行うことで作成されたファイルという感じでしょうか？（・・・浅い理解）

python 3.7に紐づいているからpython 3.6の仮想環境ではうまく動かなかったのかもしれませんね！（平然と誤情報を垂れ流すよ！!）

まとめ

以上、前回に引き続きPyMOLのCUIを利用するための方法についてでした。 適当な記事を書くとすぐに解決策が教えていただけるのでオープンソースのコミュニティに感謝感謝です。ありがとうございました！！

今回も適当調べで誤情報多そうなのでご指摘いただければ幸いです。ではでは！

*1:ソフトリンク

*2:シンボリックリンク (symbolic link)

*3:Distutilsの紹介

*4:Extensionクラス

*5:python 3 cython compiled files extension

*6:深入りしないCython入門

*7:macOS_Wikipedia