magattacaのブログ

日付以外誤報

活性化合物群を分析して指標の閾値を決めようとする話

前回、約13万個の化合物を含むEnamine_Premium_collectionから、分子量300を基準とした分子の絞り込みにより約11万個まで化合物数を減らすことができました。

まだまだ多い・・・と思っていたら@iwatobipen先生に分子を絞り込む指標を考える上でとても参考になる記事 SAR News Np.19 をご紹介いただきました。

SAR newsは日本薬学会 構造活性相関部会のHPの発行しているニュースレターで、だいたい年に２本のペースで発行されており、部会に所属していなくても過去のニュースレターを含めて自由に閲覧可能となっています（バックナンバーはこちらから閲覧可能です）。

SAR Newsで用いられている指標

さて、ご紹介いただいたニュースレターの 株式会社理論創薬研究所 吉森篤史 氏による記事では「標的タンパク質指向型フラグメントライブラリの構築」を適用例として、RDKitを用いたフィルタリングの方法が実際のコード含めて解説されています。

用いられているフィルタリング基準は、フラグメント指標「Rule of Three」で、具体的には以下となっています。

指標 基準 記事中のコード 現在のコード？
① 分子量 ≦300 AvailDescriptors.descDict['MolWt'](mol) Descriptors.MolWt(mol)
② cLogP ≦3 AvailDescriptors.descDict['MolLogP'](mol) Descriptors.MolLogP(mol)
③ 水素結合供与体数 ≦3 AvailDescriptors.descDict['NumHDonors'](mol) Descriptors.NumHDonors(mol)
④ 水素結合受容体数 ≦3 AvailDescriptors.descDict['NumHAcceptors'](mol) Descriptors.NumHAcceptors(mol)
⑤ 回転可能結合数 ≦3 AvailDescriptors.descDict['NumRotatableBonds'](mol) Descriptors.NumRotatableBonds(mol)
⑥ 極性表面積 ≦60 AvailDescriptors.descDict['TPSA'](mol) Descriptors.TPSA(mol)

記事が2010年のもののためか、現在のRDKitのモジュールを眺めてもAvailDescriptorsというものが見つけられませんでした。「化学の新しいカタチ」さんの記事「 RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ」 に同じ目的と思われるモジュール（Descriptors）が紹介されていましたので、こちらも上記の表に追加しました。

上記、指標の基準はフラグメント（"活性が弱いながらも適切な結合を形成できる小分子"）を対象としているためか、かなり小さな値となっているように見えます。今回の創薬レイドバトルの目標は、ある程度活性が強いものを取得することなので、この基準をそのまま適用するのは良くなさそうです。

また、他の指標としてリピンスキーの法則を用いる、というのも良さそうですが、前回活性化合物群の分子量のデータを眺めた印象では、分子量 500以上の大きなものも多く、すでにRule of Five を外れてしまっています。

一般的な指標をそのまま適用するだけでは、既知活性化合物群ですら外れてしまいそうなので、まずは活性化合物群の値がどのようになっているかを検証したいと思います。

PD/PD-L1阻害活性　化合物群の指標の分析

指標の計算

まずは必要なものをimportします。

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, Descriptors, PandasTools
import pandas as pd
 
print('rdkit version: ',rdBase.rdkitVersion)
print('pandas version: ', pd.__version__)
#rdkit version:  2018.09.1
#pandas version:  0.23.4

前回同様、創薬レイドバトルの共闘用シェアデータを使用します。

#csv形式のデータをPandasのDataFrameで読み込む
df = pd.read_csv('../pd1_inhibitor_dataset-master/PD1_inhibitor_dataset.csv')

#smilesカラムのSMILES情報を使ってをMolオブジェクトを追加する
PandasTools.AddMoleculeColumnToFrame(df, 'smiles')

#ROMolカラムに追加されているMolオブジェクトを使って記述子を計算し、値を追加する。
#①分子量
df['MW'] = df.ROMol.map(Descriptors.MolWt)
#②LogP
df['MolLogP'] = df.ROMol.map(Descriptors.MolLogP)
#③水素結合供与体数
df['NumHDonors'] = df.ROMol.map(Descriptors.NumHDonors)
#④水素結合受容体数
df['NumHAcceptors'] = df.ROMol.map(Descriptors.NumHAcceptors)
#⑤回転可能結合数
df['NumRotatableBonds'] = df.ROMol.map(Descriptors.NumRotatableBonds)
#⑥極性表面積
df['TPSA'] = df.ROMol.map(Descriptors.TPSA)

計算した指標の要約統計量を求めてみます。

#計算した指標のリスト
Descriptors_list = ['MW', 'MolLogP', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA']

#要約統計量を計算。小数点以下の桁数が多いので丸める
df[Descriptors_list].describe().round(2)

f:id:magattaca:20190120180844p:plain

相関係数は以下のようになっています。

#相関係数の計算
df[Descriptors_list].corr().round(2)

f:id:magattaca:20190120180932p:plain

ざっと見る限りMolLogPのみが他の指標との相関係数が負かつ絶対値が小さくなっています。

指標をプロット

数字を見てもイマイチ感触がつかめません。前回分子量を眺めた際に描画すると値の分布がわかりやすかったので今回も描画を行います。一つ一つみるのは面倒なのでSeabornを使ってペアプロットを作成します。

import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline

#Seabornを用いたペアプロットの作成
sns.pairplot(df[Descriptors_list])

f:id:magattaca:20190120181042p:plain

確かに、MolLogPだけが他の指標とのプロットした際に相関が悪そうです。

見やすさのためMolLogPを外して再度プロットして見ます。

#計算した指標のリストから MolLogPを除いたもの
Descriptors_list_2 = ['MW', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA']
sns.pairplot(df[Descriptors_list_2])

f:id:magattaca:20190120181154p:plain

残りの指標間は直線上に乗っているようにみえわかりやすくなりました。

外れ値を取り除く

各プロットを見ると、最大値となっている分子が１つだけ他の分子から遠く、外れ値のようになっています。特にNumRotatableBondsで顕著にみえるので、こちらの最大値がどのような分子なのか眺めて見たいと思います。

#NumRotatableBondsが最大値となる化合物を眺める。
#まずは最大値となる化合物の行番号を取得する
print(df['NumRotatableBonds'].idxmax())
#5と出力

#Molオブジェクトの格納されているcolumnの番号を確認
print(df.columns.get_loc('ROMol'))
#6と出力

DatafFrameの(5, 6) がみたい構造のMOlオブジェクトがある場所のようです。これを描画します。描画の手順は省きます。

f:id:magattaca:20190120181634p:plain

環状ペプチドにさらに長い鎖状構造がつながっています。他の構造と比べてかなり特異な構造なので、こちらは外れ値として省いてしまっても良さそうです。

新しいDataFrameを作成し、再度ペアプロットを作成します。

#行番号5を除去する
df_new = df.drop(5, axis=0)
sns.pairplot(df_new[Descriptors_list_2])

f:id:magattaca:20190120181824p:plain

かなり分布がわかりやすくなってきました。いずれの指標においても大きく２つのグループに分子を分けることができそうです。

これは・・・クラスタリングというやつが使えるのではないか？？？

クラスタリングによる分類

Morganフィンガープリントとタニモト係数による類似性評価

「化学の新しいカタチ」さんのこちらの記事（RDKitとk平均法による化合物の非階層型クラスタリング ）を参考に化合物を2つのクラスターに分けたいと思います。以下の操作は上記記事の処理手順をほぼそのまま使用させていただいています。

各分子のMorganフィンガープリントを計算し、分子間の類似性をタニモト係数によって比較、クラスタリングに用いるための距離行列を作成します。

from rdkit import DataStructs
from rdkit.Chem import AllChem
import numpy as np

#計算に使用するためMolオブジェクトのリストを作成
mol_list = list(df_new['ROMol'])
print(len(mol_list))
#52個の分子を含む

#Morganフィンガープリントの計算
morgan_fp = [AllChem.GetMorganFingerprintAsBitVect(x,2,2048) for x in mol_list]

#総数52の分子間で距離行列の計算
dis_matrix = [DataStructs.BulkTanimotoSimilarity(morgan_fp[i], morgan_fp[:52],returnDistance=True) for i in range(52)]

### numpy.ndarrayへの変換
dis_array = np.array(dis_matrix)
dis_array.shape

52 x 52 の距離行列が得られました。

k平均法によるクラスタリング

k平均法で２つのクラスターに分類します。

from sklearn.cluster import KMeans
import pandas as pd
kmeans = KMeans(n_clusters=2)
kmeans.fit(dis_array)
pd.value_counts(kmeans.labels_)

f:id:magattaca:20190120182815p:plain

10個の分子を含むクラス0と、42個の分子を含むクラス1に分類できました。 このクラス分類もDataFrameに追加しておきます。

class_list = []
for i in range(len(kmeans.labels_)):
    if kmeans.labels_[i] == 0:
        class_list.append('class_0')
    elif kmeans.labels_[i] == 1:
        class_list.append('class_1')
        
df_new['kmeans_class'] = class_list

ペアプロットへのクラス分類の反映

先に作成したpairplotをkmeans法で分類したクラスに応じて色分けするとどうなるか、見て見たいと思います。

もともと他の指標との相関が悪そうだったMolLogPについても、色分けしたらどうなるか気になるので、今回はpairplotに加えたいと思います。

#計算した指標のリスト（MolLogPを含む）にkmeansで形成したクラスも加える
Descriptors_class = ['MW', 'MolLogP', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA', 'kmeans_class']
sns.pairplot(df_new[Descriptors_class], hue='kmeans_class')

f:id:magattaca:20190120183026p:plain

Morganフィンガープリントをつかって２つのクラスターに分けた結果で色分けすると、指標の分布でも比較的良好に２つのクラスターに分かれているように見えます。

とくに興味深いのが、MolLogPで、ヒストグラムでは２つのクラスターが重なっていますが、ペアプロット上では割と良好に分布が分かれているように見えます。

Morganフィンガープリントがどのように計算しているのか、まだ理解していないのですが指標の分布をうまく説明するようなクラスタリングができているのは面白い結果です。

各クラスターに含まれる構造のインタラクティブな描画

プロット上での描画

各点の構造がどのようなものか、分布の広がりがありつつ２つのクラスターが綺麗に分かれているMolLogPとTPSAの２軸によるプロットを用いて可視化して見たいと思います。

@iwatobipen先生の記事 Plot Chemical space with d3js based libraryを参考にさせていただきました。

import mpld3  #別途インストールしておかないとimportできない
from mpld3 import plugins
mpld3.enable_notebook()

#描画を作成するためのモジュール
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

#svgを格納するための空のリスト
svgs_list = []


for i in range(len(df_new)):
    view = rdMolDraw2D.MolDraw2DSVG(300,300)
    
    #Molオブジェクトを含むカラムの番号
    ROMol_col = df_new.columns.get_loc('ROMol')
    
    #i行目のMolオブジェクトを取得
    mol = df_new.iloc[i, ROMol_col]
    
    #描画する分子の前処理
    prepro_mol = rdMolDraw2D.PrepareMolForDrawing(mol)
    
    #分子をコンテナに追加
    view.DrawMolecule(prepro_mol)
    
    #コンテナをファイナライズ
    view.FinishDrawing()
    
    #コンテナに書き込んだデータを取り出す
    svg = view.GetDrawingText()
    
    #データをリストに追加
    svgs_list.append(svg.replace('svg:,',''))

以上で、見たい化合物の構造情報（svg形式）を含むリストが作成できました。あとは、プロットと組み合わせて、プロット上の各点がどのような構造に相当するのかを眺めます。

fig, ax = plt.subplots()
ax.set_xlabel('MolLogP')
ax.set_ylabel('TPSA')
points = ax.scatter(df_new['MolLogP'], df_new['TPSA'])

tooltip = plugins.PointHTMLTooltip(points, svgs_list)
plugins.connect(fig, tooltip)

f:id:magattaca:20190120183609g:plain

y軸（TPSA）値が 500よりも上に位置するのがclass_1、500よりも下に位置するのがclass_0となっています。

カーソルをうまくGIFに写すことができずわかりづらくなってしまいましたが、プロット上を右上から反時計回りに順にカーソルを移動させ、構造を表示させています。 class_1は全て環状ペプチド、class_0はほとんどが環状ペプチドではない低分子ですが、左下に位置する３つは環状ペプチド構造となっています。

スライドバー形式での描画

少しわかりづらかったので、もうひとつのインタラクティブな描画方法を試します。こちらはRDKitのBlogで紹介されていたものです。

まずは　class_0 のみを要素とするDataFrameを作成し、さらにMolオブジェクトのリストとします。

df_class_0 = df_new[df_new['kmeans_class']== 'class_0']
class_0_list = list(df_class_0['ROMol'])

以下の方法ではipywidgetsを使って可視化します。

from ipywidgets import interact,fixed,IntSlider
import ipywidgets

#class_0 に分類された構造を順番に眺めていく
def class_0_viewer(idx):
    mol = class_0_list[idx]
    return(display(Draw.MolToImage(mol)))
    
interact(class_0_viewer, idx=ipywidgets.IntSlider(min=0,max=len(class_0_list)-1, step=1));

f:id:magattaca:20190120184006g:plain

class_0 のうち idx = [29, 30, 32] は環状ペプチド構造となっており、その他の構造とは少し違った形となっています。

ついでにclass_1もやってみます。同様の操作なのでコードは省略します。

f:id:magattaca:20190120184131g:plain

順番に構造をおくるとかなりわかりやすくなりました。

環状ペプチドと低分子とで２つに再度クラス分け

環状ペプチドと、それ以外で分けた方がより良いクラス分けとなりそうです。 そこでclass_0 に含まれる環状構造 をclass_1に移動し、新しいデータフレームを２つ作成します。

df_peptides = pd.concat([df_class_1, df_class_0.iloc[[29,30,32]]])
df_small_molecules = df_class_0.drop([40,41,43])  # idx=[29, 30, 32]番目の行のindexは[40,41,43]

print(len(df_peptides))  # 13と出力
print(len(df_small_molecules)) # 39と出力

合計52（13 + 39）個の２つのデータフレームに分けることができました。　処理の最初の段階で、外れ値とも思える大きな分子を一つ取り除いているので、もともとのデータセット（53個）よりも一つ少なくなっています。

活性化合物群に基づき指標の閾値を考察

さて、いろいろと操作してきたため、本来の目的を見失ってきました。

もともとの目的は、化合物ライブラリから化合物の数を絞り込むための指標の閾値を、活性があると報告されている化合物群のデータをもとに決めよう、というものでした。

活性化合物群を可視化した結果、２つのクラスターに分類することができ、これをもとに環状ペプチドと、それ以外の低分子との２つに分けることができました。このうち、「PD/PD-L1の相互作用を阻害する低分子を見つける」という目標として参照するのに適しているのは、後者となりそうです。

以下のテーブルでは、低分子化合物クラスター(small_molecules)と環状ペプチドクラスター（peptides）における各指標の最小値と最大値、およびLipinski の Rule of Five の値を並べています。

化合物を絞り込む際に用いる指標の閾値としては、値の範囲よりも少し広めにとって最右列のフィルタリング目安とするのが手かもしれません。

指標 small_molecules
min		 small_molecules
max		 peptides
min		 peptides
max		 Lipinski
Rule of 5 		フィルタリング目安
① 分子量 MW 215 835 491 2096 <= 500 > 300
（実施済）
② cLogP MolLogP -2.5 9.5 -6.35 3.20 <= 5 -3 ~ 10
③ 水素結合供与体数 NumHDonoes 1 7 9 22 <= 5 < 10
④ 水素結合受容体数 NumHAcceptors 3 9 10 26 <= 10 < 10
⑤ 回転可能結合数 NumRotatableBonds 2 15 3 39 < 15
⑥ 極性表面積 TPSA 48 250 237 713 < 250

念のため今回クラス分類した２つのデータフレームを出力しておきます。

properties = df_small_molecules.columns.values
Chem.PandasTools.WriteSDF(df_small_molecules, 'small_molecules.sdf', molColName='ROMol', idName='compound_id', properties=properties)
Chem.PandasTools.WriteSDF(df_peptides, 'peptides.sdf', molColName='ROMol', idName='compound_id', properties=properties)

まとめ

以上、本記事では分子量以外の化合物絞り込みの基準を探すべく、活性化合物群の情報を分析して見ました。

Morganフィンガープリントとk平均法を組み合わせた分類では、環状ペプチドのような構造もいくつかは鎖状低分子と同一のクラスにふくまれるなど、想定外のクラス分けがされて興味深い結果になりました。 また、ある程度構造が限られた化合物セットであれば順番に表示させるような機能を利用することで、化合物セットの雰囲気をつかむことができることがわかりました。

本記事を作成するにあたって、さまざまなページを参照させていただきました。また、多くのコードを拝借させていただいております。誤り等ございましたらご指摘いただければ幸いです。

追記

前回の記事で分子量で絞り込んだあとの化合物群（約11万個）の指標を計算してみたら、上記のテーブルのフィルタリング目安ではほとんど絞り込みにならなさそうです。 もっとちゃんと考えないとダメそうです・・・私の休日・・・

f:id:magattaca:20190121000021p:plain

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

参考記事

① SAR News No.19
Chemoinformatics Toolkits を用いた創薬システム開発におけるラピッドラタイピング　吉森　篤史
http://bukai.pharm.or.jp/bukai_kozo/SARNews/SARNews_19.pdf

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

③ 化学の新しいカタチ
RDKitとk平均法による化合物の非階層型クラスタリング

④ IS LIFE WORTH LIVING? Plot Chemical space with d3js based library #RDKit #Chemoinformatics | Is life worth living?

⑤ The RDKit Blog
RDKit: Using the new fingerprint bit rendering code

⑥ ipywidgetsのインストール方法
ipywidgetsの使い方

⑦ ipywidgets の User Guide
ipywidgets — Jupyter Widgets 7.4.2 documentation