Look Back at 2019

I’ve started the blog site from Aug. 2012. So it took over 7 years. And in this year I could get many access and comments, thanks so much. For example number of site visitors and page views are 1.5 times greater than last year.

Following KPT is just memorandum for myself .

Keep. Keep output my own activity it’s important for having fruitful discussions with any other researchers I think.

  1. I could have chances to contribute py4chemoinformatics which is a thin book for chemoinformatics and psikit which is thin wrapper for psi4. English version of py4chemoifinformatics can be found @joofio’s repository. Thank joofio for making a great contribution to Japanese <=> English translation.
  2. I could participate RDKit UGM 2019 and could get many comments about my blog post from other participants. It is a great honor for me. ;)
  3. Now I’m writing my blog post in English. It’s hard to write post in English for me but I’m going to keep it most of the posts.
  4. I started 30-60 min Running two or three times a week for my health.

Problem.

  1. Tsundoku (a Japanese term). Tsundoku means acquiring reading materials but letting them pile up in my home without reading them. In this year I bought many books about programming, statistics and printed many articles about medicinal chemistry and computational chemistry and process chemistry… I’m trying to reduce tsundoku in the winter vacation.
  2. Write code without test. I made many problems in coding because I didn’t do Test Driven Development.

Try.

  1. Improve drug discovery project process with combination of medicinal and computational chemistry.
  2. Contribute open source software development.
  3. Contribute mishima.syk, which is great community of chemo and bio informatics.
  4. Change myself. ‘It’s easy to change myself compared to change others!’

I want to wish you and your family a beautiful and happy New Year.

AI driven retrosynthesis and cooking #chmoinformatics? #cooking #pytorch

‘AI drug discovery’ is hot topics in pharma in these days. There are lots of publications about not only de novo drug design but also synthetic planning. Synthetic planning is called retro synthesis which is method to decompose target molecule to small fragments which can purchase from supplier.

For AI driven retro synthesis, input data is target molecule structure. Then the AI deradate compound to small fragments with learned chemical transformations.

BTW, do you like cooking? I like it but not good at cooking. I often search recipe with google. I think cooling and organic synthesis is same because they have recipe(experimental procedure) and it is required know how. So I wonder that are there any research about de novo cooking.

I found it !!!!

Researchers in Facebook developed inverse cooling AL with the state of the art machine learning technique.

You can find the details and code from following URL.
https://arxiv.org/abs/1812.06164
https://github.com/facebookresearch/inversecooking

The code detects photo of food and then predicts ingredient and procedure with RNN.

I felt the process is same as retro synthesis with AI.

So I had interested the code and tried to use it.

I borrowed the code from original repository and changed food photo.

The result was below.

I tried curry rice, sushi and bread. The model failed to predict curry rice and sushi.

AI missed salmon to carrot. lol

State of the art Resnet is used the algorithm but it is not perfect. It’s same as AI driven retro synthesis I think.

Current DL based model is not perfect so sometime it makes disappointing us but it has room of improvement. Near the future, we might eat food which is designed by AI. ;)

New visualization way for SAR #RDKit #chemoinformatics

Recently I found very cool visualization code for SAR that is provided by @simonduerr. It seems cool!

You can get the code from his repository.
https://github.com/duerrsimon/substrate-scope-plot

I’ve never seen such a cool visualization. Current version of the package can fix the position of the rgroup automatically.

So I tried to use the code for my own dataset.

I got SAR data from ChEMBL, its target is MPS1. At first import packages and load data.

from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import PandasTools
from rdkit.Chem import Draw
from radialscope import RadialScope as rs
from radialscope import draw_with_indeces
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
import svgutils.compose as sc
from IPython.display import SVG, HTML
import numpy as np
import pandas as pd
df = pd.read_csv('./data/mps.csv', sep=';')

Then define core and template for data extraction and visualization. Next retrieve molecules and R group information which has defined core structure.

core = Chem.MolFromSmiles('c2cn(C1CCCCC1)c3ncncc23')
templ =  Chem.MolFromSmiles('c2cn(C1CCCCC1)c3nc(C)ncc23')
matchmols = []
AlogP = []
std_val = []
for idx, mol in enumerate(mols):
    if mol.HasSubstructMatch(core):
        matchmols.append(mol)
        AlogP.append(df.AlogP[idx])
        pIC50 = np.round(9-np.log10(df['Standard Value'][idx]),2)
        std_val.append(pIC50)

rg = []
for mol in matchmols:
    sc = Chem.ReplaceCore(mol,core)
    sc = Chem.MolToSmiles(sc).replace('[1*]', '~')
    rg.append(sc)

Then, define drawing settings, I defined white_cutoff to zero because I would like to parameter with white.

settings={ 'SMILESSTRING':Chem.MolToSmiles(templ),  # use methyl groups for your rests, smiles strings can't use labels
  'use_bw_atom_theme':False,  # draw all atoms black
 'use_bold_font': True, # replace all fonts wherever possible with bold text
 'white_cutoff': 0,  #make text of labels white if greater than this number
  'scalefactor':0.7  #scale the total plot in the end, for large molecule you will need to decrease this or make the viewbox bigger in a vector software such as Inkscape
}

Next, define R group information. I used pIC50 and AlogP as SAR information. And other setting is same as original example.

Attach_atom_id is key parameter, it defines where SAR information place.

atom index 0 is dummy methyl group of the template.

radial_scope_setup_R1= {
    'rest_label':"R$_1$", # Label of the atom in the radial scope plot
    'no_wedges':len(rg), # Number of wedges in the Radial scope Plot
    'coverangle_wedges':180, # Degrees of a full circle the Rscope should cover
    'startangle':-30, #Start angle of the Rscope 
    'CMAPINNER':"Blues",  # Colormap of the inner circle. Max value has full color, options see below
    'CMAPOUTER':"RdPu", # Colormap of the outer circle. Max value has full color
    'OUTERLABEL':"pIC50", # Label of the outer circle
    'INNERLABEL':"AlogP", # Label of the inner circle
    'value_inner_circle':AlogP,
    'value_outer_circle':std_val,
    'rounding':True, # ENABLE rounding if you want something like >99 
    'rounding_boundary':10,  # cutoff for displaying > 
    'value_groups':rg, # Labels for the outer circle, you can use math using $, any ~ will be interpreted as smiles string
    'attach_atom_id': 0,
    
}

Now almost there, type following code.

scope_plot=rs(settings, radial_scope_setup_R1) # add all radial_scope dictionaries to this call
SVG('substrate_scope_replaced.svg')

It works fine!

The package can visualize multiple substituents information. However I think it should be visualize one R group information per image. Because this way is difficult to understand combination of R groups.

I don’t understand the details of the code (it’s too hacky for me) but really useful for MedChem I think.

Thanks again for sharing the nice coode.

Today’s code was uploaded my gist and github repo.

https://github.com/iwatobipen/substrate-scope-plot/blob/master/MPS1.ipynb

wet実験で大事なことってなんですかね? #souyakuAC2019

今日は12月23日!サンタさん来てくれるかなーワクワクしますよね。
ドキドキして今日はサンタさん見るまで眠れなそうなiwatobipenです。

皆さん、元気に実験していますか?私は天然物全合成バリバリのラボで年越しNMRとか年越し実験とかウキウキしながらやっていた暗い過去があります。

有機合成Wetバリバリの方と話すと意外と盛り上がるのが、やばい実験ランキングだったりしませんか。私はマスターの最後の頃、ちまちま原料合成するのがどうにも間に合わなくて、四塩化チタン1Lくらい使ってTebbe反応仕込み始めたことがあります。年末の誰もいない実験室で一人焦って雑に作業した結果、シリンジが詰まって半泣きなったことや、金属ナトリウム50g位を使ってBarch還元したときに反応が暴走仕掛けた苦い記憶があります。他にも大量のLiAlH4還元も、、、今となっては絶対やりたくないですね。

さて、長い前置きでしたが、あなたにとってWet実験で大事なことってなんですか。再現性?コスト? いえ一番優先されるべきは安全性ですよ! 実験には危険がいっぱい。安全を甘く見ていると後で手痛いしっぺ返しを喰らいます。自分の健康を損なうことも。企業であれば操業停止などもありえます。

有機合成で最近ニュースになった大きな話題は2008年に起こったSangjiさんの事故かと思います。t-BuLiを扱う際に起きたこの痛ましい事故は改めて実験時の安全管理大切さを考えるきっかけとなりました。
https://cen.acs.org/articles/87/web/2009/01/Researcher-Dies-Lab-Fire.html

なぜこんな記事を書いているかというと、最近気になる文献を読んだからです
https://www.nature.com/articles/s41557-019-0375-x

先の事故以降あまり、死亡に至る事故の記事を目にしてこなかったので驚いたのですが、この文献のTable1に死亡事故まで至ったアクシデントのリストがあり、2008年の事故以降も事故は起きています。有毒ガスや、引火性の期待に起因する事故が多いように見えます。ガスは目に見えませんし、知らぬ間にドラフトから漏洩するリスクもありますね。

このような事故が起きる要因として著者らは事件の安全性における包括的なデータセット、それに基づく教育の仕組みがないことが要因だろうと解析しています。たしかに私も大学の講義でLab safetyナド教わったことはありませんね、、、。水素添加時の発火とか、やべー火い吹いたわ、てへっ。みたいなヒヤリで済むとその場限りなっちゃいます。こういったヒヤリの事例はかなりの頻度で起きているのではないでしょうか。それを元に再発を防ぐことが大事ですよね。

ですが、そのような事例をまとめて教育する仕組みって意外とないですし、教育をどれだけしても完璧とは言えません。なぜなら多くの実験は人がするものでありそれが故に、リスクがあるのです。次の一文がそのことを指摘しております。

“Of all the variables in accident prevention, the human behaviour
variable, even with education, was the hardest to control”

更に事故が起きるリスクは若手の研究者のほうが高い傾向があります。これはやはり経験や知識がまだ未熟であるためでしょう。ただし、経験が多いから安全ということはなく慢心は思わぬ事故を招きます。常にリスクを評価し実験に望むことが大事です。

さて、安全が大事、安全第一というのはわかっているのに、意外と軽視される理由は何なんでしょう?著者らはBarriers to safety research というセッションでそのことについて議論しています。読んでいてそうだよなーって思ったのは、下記のパートです

The issue of ‘academic freedom’ is often raised as an objection to safety practices. One study found that 15% of researchers believed that safety regulations interfered with productivity, and 23% believed that they impeded the scientific discovery process.

安全対策とか教育とか、保護具とかあれこれ対策するのは時間もコストもかかります。世界初の結果を出して論文投稿して、サイエンスを切り開くことこそが研究者の生きる道。安全は二の次!となっちゃうとよろしくない。面倒なのはわかります。でもそれでやらなくていいことにはならないんですよね。

このアーティクルはあくまでアカデミアでのLab safetyという切り口で議論していますが、企業の研究者の皆さんどうですか?

色々自分の身を守るための安全対策、マニュアル、その他もろもろ、煩わしいなって思ったりしませんか。

研究者に必要なのは安全顧みず突き進む度胸ではないのです。猪突猛進に研究にのめり込むのは大切。でも安全あってこそですよ。怪我してからでは遅いのです。

幸い私は大きな怪我もなくここまで研究者生活を遅れています。それは、安全にすごい気を使っていたからではなく、多分たまたまなんです。労働災害に関する有名な法則にハインリッヒの法則があります。ヒヤリで済むか、重大事故になるかは意外と紙一重かもしれません。

最近けっこう面白い論文見つけてそれについても書こうかなと思っていたんですが、だいぶ趣の変わった文章を書いてしまいました。

今年も終わりに近づいてきました。今年もこれからもずっと皆様が安全に、素敵な研究ライフを遅れることを心より願い今日の駄文を終わりにしたいと思います。

最後まで読んでくださりありがとうございました☆。

Clustering molecules with HDBSCAN and predict chemical space with scikit-learn #chemoinformatics #RDKit

The contents of the post is almost same as yesterday’s one that was for souyaku-advent calendar 2019.

In drug discovery area we often discuss about compound chemical space. Chemist should be explore chemical space to find new IP space etc. It is defined with descriptors, fingerprint and biological response etc… It is often high dimensional space so we would like to reduce the dimension with some mathematical approach such as PCA, TSNE and UMAP.

And also we use clustering method for grouping molecules.

There are many clustering method for example k-means, hierarchical clustering. However k-means need to define k value and hierarchical method need to define threshold for doing that. So try and error is required. If you would like to get clustering label conveniently, I recommend to try to use HDBSCAN. HDBSCAN clusters data with their density. Ref URL is below.
https://link.springer.com/chapter/10.1007%2F978-3-642-37456-2_14

Fortunately we can use HDBSCAN from python. It is easy to install via pip or conda. Please check original document to install. ;)

I tried to demo with HDBSCAN. Following data is bollowed from chaper8 in py4chemoinformatics (Thanks@fmkz___ for providing data and code!)

Let’s write code! Import packages at first. I’ll use mlinsght later.

%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

I defined tanimoto-dis tfunction for HDBSCAN. HDBSCAN requires metric arg. for measuring distances of each data point. And there isn’t tanimoto distance as a default metric. If user would like to original metric function, pass ‘pyfunc’ for metric arugment and your own function is passed to func argument. Func should be took two numpy array and return distance.
@@@@@
My tanimoto dist function is tooooooooooooooooooooooooo slow it has room for improvement I fill improve it later. ;)
@@@@@

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)

Now ready, let’s check data without and with HDBSCAN. I made two plots one is coloured with company name and the other is coloured with HDBCAN label.

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)

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)
coloured with company name
coloured with hdbscan label.

Coloured with Gray color point is not labelled data. The result shows HDBSCAN could do almost reasonable labelling in this case.

Plot new compound with defined chemical space

Next topic is how to plot new compound in pre-defined chemical space.
Theoretically, PCA and tSNE can’t add new data to pre-defined latent space. Re calculation is required if user would like to add new data. But we would like to map new designed compounds to current chemical space. Hmm…. how to do it?

The answer for the question is provided from mlinghts. The package has PredictableTSNE class. The class take transformer and estimator object. And it makes predictive model for transformer output (such as PCA, TSNE output) with estimator. It will be clear if you read source code. After learning, PredictableTSNE can predict latent space from data features.

Following code, I used random sampling and random forest and GP regressor for estimator.

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))])
blue is training data, orange is test data.

Use predictiveTSNE for new compound chemical space prediction. Following examples are in RF case and GP case.


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)


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))
RF as an estimator
GP as an estimator

In above case, random forest seems better performance for new compound prediction.

The method of estimator has a big influence for preditiveTSNE performance. It is useful for plot new data but it has risk for miss-leading for design.

At the end, I checked cluster label and preditiveTSNE data.

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


x training data, o test data, colour : class label.

Not only training data but also test data which has same colour are plotted almost same area. It indicates that the approach works well in this case.

Finally we often discuss chemicals pace and use the term ‘chemical space’ but it is very difficult to define, use and compare it.

All code is uploaded in my gist. Any comments and/or suggestions are greatly appreciated.

Chemicalspace x clustering = ? #souyakuAC2019 #chemoinformatics

クリスマス、年末年始が近づいてきた今日この頃。皆様いかがお過ごしでしょうか。寒さにはめっぽう弱い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の性能にだいぶ作用されるので実戦投入時は要注意ですね。

さt,最後に比較的良さげだった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')


xはトレーニング、◯はテストデータです。色がクラスターです同じクラスターは同じあたりに投影されているのでそこそこ上手くワークしているようです。

今回完全にランダムサンプリングしたのでうまく行っていますが、全く道のデータが来ると性能はガクッと落ちるかもしれません。

俺のデザインする化合物はお前らのモデルで定義した空間に収まるような代物じゃねえぜ!という熱い思いの化合物を作らないとですねw

冗長ですが、まあ近傍のデザインをする上では面白いかなと思い紹介しました。HDBSCANやTSNEはデータが増えると計算コストがけっこうかさむので注意しましょう。

はあー出かける前に書き終わったのでホッとした。

コード全体は以下のGISTからどぞー

TIPS of rdkit #RDKit #memorandum #chemoinformatics

As you know, rdkit is very attractive and active project for chemoinformatics. Recent version of rdkit has heteroshuffle enumerator. It is useful for generate new molecules. You can find the details of the method below.
https://www.rdkit.org/docs/source/rdkit.Chem.EnumerateHeterocycles.html

However in the drug discovery project, it is not needed enumerate all aromatic rings in the molecule I think. So enumeration of targeted ring is useful I think.

Fortunately EnumerateHeterocyles uses rdkit Reaction method by using defined reaction rules. It means that specific atoms can protect with setting atom property. Let’s write code! Following code I use omeprazole as example. And I got ring information to protect specific atoms.

from rdkit import rdBase
from rdkit import Chem
from rdkit.Chem import EnumerateHeterocycles
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import AllChem
import copy
omeprazole = Chem.MolFromSmiles('CC1=CN=C(C(=C1OC)C)CS(=O)C2=NC3=C(N2)C=C(C=C3)OC')
ringinfo = omeprazole.GetRingInfo()
rings = ringinfo.AtomRings()
print(rings[0])

At first enumerate heterocycles without any restriction.

res = EnumerateHeterocycles.EnumerateHeterocycles(omeprazole)
res = [m for m in res]
len(res)
> 384

Check generated molecules.

Draw.MolsToGridImage(res[:20], molsPerRow=3)

Next, protect pyridine ring and generate heterocycle of bicyclic moiety. As expected, number of generated molecules are reduced.

protected_omeprazole = copy.deepcopy(omeprazole)
for atmidx in rings[0]:
    atom = protected_omeprazole.GetAtomWithIdx(atmidx)
    atom.SetProp('_protected', '1')
res2 = EnumerateHeterocycles.EnumerateHeterocycles(protected_omeprazole)
res2 = [m for m in res2]
len(res2)
> 96

Check the structures.

Draw.MolsToGridImage(res2[:50], molsPerRow=5)

Protect works fine!

And I checked whether the function keep 3D structure of the molecule.

Also I worked fine.

So it is useful for generate specific focused library.

I uploaded the code to gist.