クリスマス、年末年始が近づいてきた今日この頃。皆様いかがお過ごしでしょうか。寒さにはめっぽう弱いiwatobipendです。これは2019年の創薬アドベントカレンダー（DRY)の記事です。
https://adventar.org/calendars/4106
化合物空間（ケミカルスペース）は広大な宇宙だ！新しい化合物空間を求めて旅をするのがケミストだ！とか、私達のライブラリは特定の化合物空間を専有しています。とか、化合物空間という言葉良く聞きますよね。化合物空間というのは、フィンガープリントとか記述子をもとに定義される化合物が取り得る空間のことを指します。
一般にフィンガープリントなどは高次元（1024次元とか）で私のような、凡人には理解できないので、上の高次元情報を適当に圧縮して2次元とかで眺めます。うまく圧縮できると、ケミストが見た構造クラスタと出力結果がいい感じにマッチしてスキッとするわけです。更にここに活性はパラメータを入れて今度はこの辺を攻めるか！という創薬という大海原への地図ができます（俺は創薬王になる！ドンっ）
今日はこのあたりで使えそうなものをちょっと紹介したいなと思います。
いい感じクラスタリングしてほしいのよね
化合物を似ているものどうしてクラスタリングしておくと色々と便利です。化合物の傾向を見たりするときなど。K近傍法のように予め分けたいクラスタ数を指定する方法と、階層別クラスタリングのように、あとで人が閾値を決めて分ける方法ナド色々方法がありますが、もういちいちそんなんヤッテランネーヨ！ほれ最近流行りのエーアイとかなんかよろしくやれるもんないのかよ。 という方もいるかと思います（僕だけ？）。そんな方におすすめなのがHDBSCAN。本手法は階層別クラスタリングと合わせて密度が濃い部分をまとめてクラスタにするといったことを自動でやってくれます。詳しいアルゴリズムは文献を読んでください。
それでは早速手を動かしてみましょう。データはPy4chemoinformaticsで@fmkz___さんが提供してくれているものを利用します。まず必要なパッケージを読み込みます。HDBSCANはAnaconda経由で、MLINSGIHTSはPIPで入れます。後者は後ほど使います。
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Draw
from rdkit import RDLogger
from sklearn.model_selection import train_test_split
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.manifold import TSNE
from mlinsights.mlmodel import PredictableTSNE
from hdbscan import HDBSCAN
sns.set_context('poster')
sns.set_style('white')
sns.set_color_codes()
plot_kwds = {'alpha' : 0.5, 's' : 80, 'linewidths':0}
RDLogger.DisableLog('rdApp.*')
seed = 794
データの読み込みは8章ほぼ丸コピです。tanimoto_distの定義からHDBSCANのオブジェクト作りに入ります。HDBSCANは密度計算のため距離METRICSを図る必要がありますが、Chemo用パッケージじゃないのでユークリッド距離、コサイン距離などはありますが谷本距離はありません。このような場合HDBSCANのmetricの引数に’pyfunc’を指定し、使う距離測定関数をfuncで渡してあげると独自の関数で評価してくれます。なお、Funcに渡す関数は2つのNumpy arrayをもらい距離の数値を返す必要があるのでRDKITのBitVecをちょくで渡すとエラーになります。
oxrs = [("CHEMBL3098111", "Merck" ), ("CHEMBL3867477", "Merck" ), ("CHEMBL2380240", "Rottapharm" ),
("CHEMBL3352684", "Merck" ), ("CHEMBL3769367", "Merck" ), ("CHEMBL3526050", "Actelion" ),
("CHEMBL3112474", "Actelion" ), ("CHEMBL3739366", "Heptares" ), ("CHEMBL3739395", "Actelion" ),
("CHEMBL3351489", "Eisai" )]
fps = []
docs = []
companies = []
mol_list = []
for cid, company in oxrs:
sdf_file = os.path.join("data", cid + ".sdf")
mols = Chem.SDMolSupplier(sdf_file)
for mol in mols:
if mol is not None:
mol_list.append(mol)
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
docs.append(cid)
companies.append(company)
fps.append(arr)
fps = np.array(fps)
companies = np.array(companies)def tanimoto_dist(ar1, ar2):
a = np.dot(ar1, ar2)
b = ar1 + ar2 - ar1*ar2
return 1 - a/np.sum(b)
clusterer = HDBSCAN(algorithm='best', min_samples=5, metric='pyfunc', func=tanimoto_dist)
clusterer.fit(fps)
docs = np.array(docs)
trainIDX, testIDX = train_test_split(range(len(fps)), random_state=seed)
準備ができたので、次元圧縮と、クラスタリングの結果を見ましょう。今回次元圧縮にはTSNEを使います。Scikit-learnは便利ですね。
tsne = TSNE(random_state=seed) res = tsne.fit_transform(fps)
まずは会社ごとに色つけしてみてみます。きれいに分かれていますね。
plt.clf() plt.figure(figsize=(12, 6)) sns.scatterplot(res[:,0], res[:,1], hue=companies, **plot_kwds)
次にHDBSCANでクラスタリングした結果を反映させてみます。ラベルはclusterer.labels_というプロパティに格納されています。
plt.clf()
plt.figure(figsize=(12, 6))
palette = sns.color_palette()
cluster_colors = [sns.desaturate(palette[col], sat)
if col >= 0 else (0.5, 0.5, 0.5) for col, sat in zip(clusterer.labels_, clusterer.probabilities_)]
plt.scatter(res[:,0], res[:,1], c=cluster_colors, **plot_kwds)
グレーの部分はクラスターに割り振られていない部分です。Min clustersize, nodesizeなどのオプションで変わるので思った結果にならないときは適当に調整してみましょう。下の例ではなんとなく点が集まっているところが一つのクラスタになってるのでいい感じじゃないでしょうか。
新しい化合物をプロットしたいんだけどさ。。。
さて、PCAやTSNEでケミカルスペースを表現できることはわかりました。しかしながらPCAやTSNEという手法は手元のデータに対して実施して跡からデータを追加したいときは、ゼロから計算し直す必要があります。なぜなら化合物の分散や距離の比較は全化合物に対してやりたいからです。一方通常の創薬プロジェクトは常に新しく化合物が生み出されます。それらがどこにいるか知りたいわけです。SOFにもTSNEで新しいデータいい感じに入れたいんだけどって質問がたまにあります。
このニーズに対する答えの一つがmlinsightsです。このパッケージにはPredictableTSNEというクラスがあります。ソースコードを見ればわかるのですがトレーニングセットでTSNEを実施ある特徴量を圧縮した次元に射影。別途渡すEstimatorクラスで特徴量＝＞圧縮された次元の予測モデルを作ります。
これで新しいデータも圧縮後の次元の予測値が得られるわけです。渡すEstimatorはSKLEARNのものは何でもオッケーです。今回はRandom forestとGPを使ってみます。まずデータをランダム分割します。下図のようにまんべんなくサンプリングしています。
trainFP = [fps[i] for i in trainIDX] train_mol = [mol_list[i] for i in trainIDX] testFP = [fps[i] for i in testIDX] test_mol = [mol_list[i] for i in testIDX] allFP = trainFP + testFP tsne_ref = TSNE(random_state=seed) res = tsne_ref.fit_transform(allFP) plt.clf() plt.figure(figsize=(12, 6)) sns.scatterplot(res[:,0], res[:,1], hue=['train' for i in range(len(trainFP))] + ['test' for i in range(len(testFP))])
ではトレーニングデータでケミカルスペースを学習、テストデータのケミカルスペースを予測、両方プロットの順でおってみましょう。PredictableTSNEのTransformer, estimatorに任意のオブジェクトを渡せば良いです。TSNEはYを与えても無視しますが、ないとエラーになるのでダミーを入れています。モデルができたら重ねてプロットしてみます。赤がテストデータです。
rfr = RandomForestRegressor(random_state=seed) tsne1 = TSNE(random_state=seed) pred_tsne_rfr = PredictableTSNE(transformer=tsne1, estimator=rfr, keep_tsne_outputs=True) pred_tsne_rfr.fit(trainFP, list(range(len(trainFP)))) pred1 = pred_tsne_rfr.transform(testFP) plt.clf() plt.figure(figsize=(12, 6)) plt.scatter(pred_tsne_rfr.tsne_outputs_[:,0], pred_tsne_rfr.tsne_outputs_[:,1], c='blue', alpha=0.5) plt.scatter(pred1[:,0], pred1[:,1], c='red', alpha=0.5)
構造が見えないのでなんとも言えませんがまあそれっぽくテストデータも投影しているかな？
続いてGPモデルで。
gbr = GaussianProcessRegressor(random_state=seed) tsne2 = TSNE(random_state=seed) pred_tsne_gbr = PredictableTSNE(transformer=tsne2, estimator=gbr, keep_tsne_outputs=True) pred_tsne_gbr.fit(trainFP, list(range(len(trainFP)))pred2 = pred_tsne_gbr.transform(testFP) plt.clf() plt.figure(figsize=(12, 6)) plt.scatter(pred_tsne_gbr.tsne_outputs_[:,0], pred_tsne_gbr.tsne_outputs_[:,1], c='blue', alpha=0.5) plt.scatter(pred2[:,0], pred2[:,1], c='red', alpha=0.5))
あらら、、だめだめ。Estimatorの性能にだいぶ作用されるので実戦投入時は要注意ですね。
さｔ，最後に比較的良さげだったRFのモデルのデータにクラスタ情報を足してみましょう
allmol = train_mol + test_mol
fps2 = []
clusterer2 = HDBSCAN(algorithm='best', min_samples=5, metric='pyfunc', func=tanimoto_dist)
for mol in allmol:
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
fps2.append(arr)
fps2 = np.array(fps2)
clusterer2.fit(fps2)
plt.clf()
plt.figure(figsize=(12, 6))
plt.scatter(pred_tsne_rfr.tsne_outputs_[:,0], pred_tsne_rfr.tsne_outputs_[:,1], c=clusterer2.labels_[:len(trainFP)], alpha=0.5, marker='x')
plt.scatter(pred1[:,0], pred1[:,1], c=clusterer2.labels_[len(trainFP):], alpha=0.5, marker='o')
ｘはトレーニング、◯はテストデータです。色がクラスターです同じクラスターは同じあたりに投影されているのでそこそこ上手くワークしているようです。
今回完全にランダムサンプリングしたのでうまく行っていますが、全く道のデータが来ると性能はガクッと落ちるかもしれません。
俺のデザインする化合物はお前らのモデルで定義した空間に収まるような代物じゃねえぜ！という熱い思いの化合物を作らないとですねｗ
冗長ですが、まあ近傍のデザインをする上では面白いかなと思い紹介しました。HDBSCANやTSNEはデータが増えると計算コストがけっこうかさむので注意しましょう。
はあー出かける前に書き終わったのでホッとした。
コード全体は以下のGISTからどぞー
Is life worth living? 