化合物生成をディープラーニングで？

![banner.png](https://qiita-image-store.s3.amazonaws.com/0/259500/d85e0519-24cb-d85a-4622-fc7a51b3c2cf.png)

#はじめに
過去に書いたものを再掲しています。
雑記に近いかもしれません。あしからず。

#新規分子生成？
新しい分子を生成する。特に「欲しい物性を持つ有用な新規分子」の設計について、今までどうしてきたのでしょうか？例えば、創薬の分野では化学の基礎理論は勿論のこと、経験則、Tanimoto係数 、量子化学計算　etc....（他にもあると思いますが、）を用いて創出していると思われます（自分調べ）。


最近、機械学習に上記のことを行わせようという大きな動きがあるようです。


その中で今回注目したのが、引用数が213と多く、ハーバード大のAspuru-Guzikグループが発表した
「Automatic Chemical Design Using a Data-Driven Continuous Representation of Molecules」(1)
という論文とこれを元に作られたプログラム
「Chemical VAE」(2)です。
![banner.png](https://qiita-image-store.s3.amazonaws.com/0/259500/d85e0519-24cb-d85a-4622-fc7a51b3c2cf.png)

#Chemical VAEって何？

これは、SMILES2vecというword2vec（Seq2Seq）を利用した技術です。これについては前回の[私の記事](https://qiita.com/muuu4649/items/7dcbe68db453b1311a6b　"私の記事")を見ていただければ幸いです。

以下は論文を読んで、自分なりに考えた新規分子生成の流れです。

まず、化合物の表す文字列をエンコーダーでベクトル化し、潜在空間（ベクトル空間）を生成させます。このベクトル空間の位置１点１点はSMILESの文字列であり、位置（後にzと表されます。）が近いほど構造が近いものが存在するそうです。それをデコーダーでできるだけ、同じような文字列に戻します。また、エンコードとデコードがうまく行くように、エンコーダーとデコーターを学習をさせていきます。


![smv.jpg](https://qiita-image-store.s3.amazonaws.com/0/259500/d24ad604-1eea-e6bd-c857-cb88bf665a8a.jpeg)

その後、（「同時に」と言ってもいいかもしれません。）ベクトル空間と対応する分子の物性値をニューラルネットワークで学習することで、下図のようにf(z)を生成します。

![スクリーンショット 2018-12-31 17.15.12.png](https://qiita-image-store.s3.amazonaws.com/0/259500/9cc0a34f-6ff0-bb02-8f1d-1e1235f9f0a6.png)


#（注）ここからは特に読解が難しかったため、予想に近いものになります。

新規分子と言っても欲しい物性がもたなければ意味がありませんよね？
今回紹介した技術は新規分子生成ではあるものの、
「物性の良い既知の分子に近い構造の分子は同様に物性が良い場合がある（多い？）」という薬学的経験則？のような考えがあるようです。
つまり、この学習した潜在空間中に物性の良い既知の分子をエンコードして、その分子の潜在空間中の位置周辺を探して新規分子を探すという過程をとるようです。その後、学習させたデコーダでその分子のSMILES生成します。最後に生成する際、分子として成り立つかどうかをRDkitで判別させます。

![oc-2017-00572f_0005.gif](https://qiita-image-store.s3.amazonaws.com/0/259500/48602be8-92f9-b866-ca8e-44d53562a55e.gif)

下図はその結果のようです。四角で囲んでいるものが中心分子です。そこから潜在空間の位置関係が表されています。

![Figure2.jpg](https://qiita-image-store.s3.amazonaws.com/0/259500/8107117c-b3a9-fd3f-60f6-39465601b1ee.jpeg)

ここまでの流れが私の理解内です。

#実際にGitHubのExampleをやってみた。
Exampleをやる際は、chemical_vaeをcondaやpipで入れてください。

まずは、import部分

```intro_to_chemvae.ipynb
# tensorflow backend
from os import environ
environ['KERAS_BACKEND'] = 'tensorflow'
# vae stuff
from chemvae.vae_utils import VAEUtils
from chemvae import mol_utils as mu
# import scientific py
import numpy as np
import pandas as pd
# rdkit stuff
from rdkit.Chem import AllChem as Chem
from rdkit.Chem import PandasTools
# plotting stuff
import matplotlib.pyplot as plt
import matplotlib as mpl
from IPython.display import SVG, display
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
```

データセットはzincのデータセットを使います。このデータセットは、SMILESと物性（QED（薬らしさ評価）,SAS(synthetic accessibility score),logP（オクタノール係数））の数値が入っています。

また、
・smiles_1は中心分子を指定します。
・noiseは潜在空間中の中心分子からの距離（z）です。
・z範囲内でランダムサンプリングをしており、SMILESが成り立つものが試行１回だけでは、みつからない場合があるため、500回for文で回してみました。
・Reconstructionは中心分子を学習したエンコーダーでベクトル化したものを学習したデコーダーで出力したものです。（下の結果はうまくいってなさそうですが、学習の仕方、パラメーターを変えるべきでしょう。）



```
vae = VAEUtils(directory='../models/zinc_properties')
smiles_1 = mu.canon_smiles('CSCC(=O)NNC(=O)c1c(C)oc(C)c1C')

for i in range(500):
   X_1 = vae.smiles_to_hot(smiles_1,canonize_smiles=True)
   z_1 = vae.encode(X_1)
   X_r= vae.decode(z_1)

   print('{:20s} : {}'.format('Input',smiles_1))
   print('{:20s} : {}'.format('Reconstruction',vae.hot_to_smiles(X_r,strip=True)[0]))

   print('{:20s} : {} with norm {:.3f}'.format('Z representation',z_1.shape, np.linalg.norm(z_1)))


   print('Properties (qed,SAS,logP):')
   y_1 = vae.predict_prop_Z(z_1)[0]
   print(y_1)
　　noise=3.0
   print('Searching molecules randomly sampled from {:.2f} std (z-distance) from the point'.format(noise))
```

・出力結果


```
Using TensorFlow backend.
Standarization: estimating mu and std values ...done!
Input                : CSCC(=O)NNC(=O)c1c(C)oc(C)c1C
Reconstruction       : CH1nCNc1Cs)Nccccc(CCc1)c3
Z representation     : (1, 196) with norm 9.901
Properties (qed,SAS,logP):
[0.72396696 2.1183593  2.1463375 ]
Searching molecules randomly sampled from 3.00 std (z-distance) from the point
```


最後は、見つけてきた、ユニークかつRDkitに判別されたものを取得してきます。


```
   
   df = vae.z_to_smiles( z_1,decode_attempts=100,noise_norm=noise)
   print('Found {:d} unique mols, out of {:d}'.format(len(set(df['smiles'])),sum(df['count'])))
   print('SMILES\n',df.smiles)
   if sum(df['count']) !=0:
      df1=pd.DataFrame(df.smiles)
   df1.to_csv("result1.csv",mode='a',index=False,header=False)

```

以下出力結果

```result1.csv
ON cCO=COCC(O)ccN2cs2c
CCCCCCNc-1cO-SCOCCcccc1
CC1CCcC(-nOcc1ccccCCC1)c1 O
OC (C)C(=Occc3cccccccc)CB
CCC1oNCc2cCcccccc2cccc1 1 1
CO CC(1c(O=O1O(1cO)nC))1
C=C1nn(=O)SnNccccccocc1C
CC C@Cs(=CN=11cccc2cc1Cc)2c1
C OcCCc(CO)c1nccc=Occccc1C
O1Cc(c1CCO)CNCC=BBOCCCN
CC ON(FCNN(C)ccc(Ocn1)1)l
C1ccnccnccccccccncccscc1
CC(CScc(c1cOn1nc1CCl)C)1
CCCCc(-ncccc21nc1c1c2)1CCC
C cnc(Cnncncc(C())Cl)Cl1 1
```

RDkitを通したにも関わらず、おそらく分子にならないものもちらほら。。。
しかし、多くのSyntax errorを回避しているようです。RDkit有能。。

以上

#参考
１）Automatic Chemical Design Using a Data-Driven Continuous Representation of Molecules
https://pubs.acs.org/doi/abs/10.1021/acscentsci.7b00572
２）chemical_vae
https://github.com/aspuru-guzik-group/chemical_vae
３）DeepLearningによる化合物生成(薬、有機発光分子)
https://ritsuan.com/blog/8480/