Smiles2vecで物性予測をしよう。

Smiles2vecとは？

簡単に言うと自然言語処理（NLP）の分野の技術で、文字列をベクトルに変換するというものです。文字列で物性予測って何?という方も多いでしょう。

Smilesとは、SMILES記法（スマイルスきほう、英: Simplified molecular input line entry system）とは、分子の化学構造をASCII符号の英数字で文字列化した構造の曖昧性の無い表記方法である。[wikipediaより](https://ja.wikipedia.org/wiki/SMILES記法　"Wikipedia")

化学構造の情報を持った文字列から原子情報だけではなく、時系列情報もベクトルに変換し、その特徴量を入力変数として、目的の物性を予測するという感じでしょうか？
それでは、やっていきます。

Smiles2vecの構造について

上でも話したようにSMILES2vec はNLPの分野の技術で、Seq2Seqという名前の技術が元になっています。この技術はRNN（リカレントネットワーク：再帰的ニューラルネットワーク）の考えに基づいたLSTM（Long short-term memory）やGRU（Gated recurrent unit）などの層を使います。LSTMやGRUはそれぞれ、時系列データのベクトル化を行うときに用いられるニューラルネットワークです。

実際にやってみた。

下記リンクを参考にさせていただきました。
https://github.com/Abdulk084/Smiles2vec
今回は　Tetrahymena pyriformis IGC50 48h toxicityについての活性を予測します。

・importとデータのセット

```smiles2vec.py
import keras
from keras import backend as K
from keras.callbacks import ModelCheckpoint
from keras import optimizers
from keras.layers import Dense, Dropout
from keras.models import Sequential
from keras.layers.recurrent import LSTM
from sklearn import datasets
from sklearn import metrics
import numpy as np
import pandas as pd
from keras import layers
from keras.layers import Input, Add, Dense, Activation, ZeroPadding2D, BatchNormalization, Flatten, Conv2D, AveragePooling2D, MaxPooling2D, GlobalMaxPooling2D
from keras.models import Model, load_model
from keras.utils import layer_utils
import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem
from keras.models import Sequential, Model
from keras.layers import Conv2D, MaxPooling2D, Input, GlobalMaxPooling2D
from keras.layers.core import Dense, Dropout, Activation, Flatten
from keras.optimizers import Adam
from keras.callbacks import ReduceLROnPlateau


data = pd.read_excel(r'IGC50.xlsx')

X_train_smiles = np.array(list(data["smiles"][data["split"]==1]))
X_test_smiles = np.array(list(data["smiles"][data["split"]==0]))
print(X_train_smiles.shape)
print(X_test_smiles.shape)
```

・SmilesをOne-hot表現に変換

```
assay = "Activity"
Y_train = data[assay][data["split"]==1].values.reshape(-1,1)
Y_test = data[assay][data["split"]==0].values.reshape(-1,1)

charset = set("".join(list(data.smiles))+"!E")
char_to_int = dict((c,i) for i,c in enumerate(charset))
int_to_char = dict((i,c) for i,c in enumerate(charset))
embed = max([len(smile) for smile in data.smiles]) + 5
print (str(charset))
print(len(charset), embed)

def vectorize(smiles):
       one_hot =  np.zeros((smiles.shape[0], embed , len(charset)),dtype=np.int8)
       for i,smile in enumerate(smiles):
           #encode the startchar
           one_hot[i,0,char_to_int["!"]] = 1
           #encode the rest of the chars
           for j,c in enumerate(smile):
               one_hot[i,j+1,char_to_int[c]] = 1
           #Encode endchar
           one_hot[i,len(smile)+1:,char_to_int["E"]] = 1
       #Return two, one for input and the other for output
       return one_hot[:,0:-1,:], one_hot[:,1:,:]


X_train, _ = vectorize(X_train_smiles)
X_test, _ = vectorize(X_test_smiles)
X_train[8].shape
vocab_size=len(charset)
```

・Smiles2vec層の定義

```
from keras.preprocessing.text import one_hot
from keras.preprocessing.sequence import pad_sequences
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Flatten
from keras.layers.embeddings import Embedding

model = Sequential()
model.add(Embedding(vocab_size, 50, input_length=embed-1))
model.add(keras.layers.Conv1D(192,10,activation='relu'))
model.add(BatchNormalization())
model.add(LSTM(224,return_sequences=True))
model.add(LSTM(384,return_sequences=True))
model.add(Flatten())
model.add(Dense(100, activation='relu'))
model.add(Dropout(0.4))
model.add(Dense(1, activation='linear'))

```

・評価とcsvへの保存

```

def coeff_determination(y_true, y_pred):
   print(y_pred)
   from keras import backend as K
   SS_res =  K.sum(K.square( y_true-y_pred ))
   SS_tot = K.sum(K.square( y_true - K.mean(y_true) ) )
   return ( 1 - SS_res/(SS_tot + K.epsilon()) )
def get_lr_metric(optimizer):
   def lr(y_true, y_pred):
       return optimizer.lr
   return lr

optimizer = Adam(lr=0.00025)
lr_metric = get_lr_metric(optimizer)
model.compile(loss="mse", optimizer=optimizer, metrics=[coeff_determination, lr_metric])

callbacks_list = [
   ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=5, min_lr=1e-15, verbose=1, mode='auto',cooldown=0),
   ModelCheckpoint(filepath="weights.best.hdf5", monitor='val_loss', save_best_only=True, verbose=1, mode='auto')

]


history =model.fit(x=np.argmax(X_train, axis=2), y=Y_train,
                             batch_size=128,
                             epochs=150,
                             validation_data=(np.argmax(X_test, axis=2),Y_test),
                             callbacks=callbacks_list)

y_pred=model.predict(np.argmax(X_test,axis=2))
df=pd.DataFrame(y_pred)
df.to_csv("result.csv",index=False,header=False)

```

・今回使用したデータセット（一部抜粋）　IGC50.xlsx

#結果

縦軸：実測値　横軸：予測値

```
Epoch 00101: ReduceLROnPlateau reducing learning rate to 3.051757957450718e-08.

Epoch 00101: val_loss did not improve from 0.27625
Epoch 102/150
1434/1434 [==============================] - 16s 11ms/step - loss: 0.2847 - coeff_determination: 0.7465 - lr: 3.0518e-08 - val_loss: 0.2815 - val_coeff_determination: 0.7069 - val_lr: 3.0518e-08

Epoch 00102: val_loss did not improve from 0.27625
Epoch 103/150
1434/1434 [==============================] - 16s 11ms/step - loss: 0.2527 - coeff_determination: 0.7762 - lr: 3.0518e-08 - val_loss: 0.2815 - val_coeff_determination: 0.7069 - val_lr: 3.0518e-08

Epoch 00103: val_loss did not improve from 0.27625
Epoch 104/150
1434/1434 [==============================] - 16s 11ms/step - loss: 0.2748 - coeff_determination: 0.7527 - lr: 3.0518e-08 - val_loss: 0.2815 - val_coeff_determination: 0.7069 - val_lr: 3.0518e-08

Epoch 00104: val_loss did not improve from 0.27625
Epoch 105/150
1434/1434 [==============================] - 16s 11ms/step - loss: 0.2793 - coeff_determination: 0.7503 - lr: 3.0518e-08 - val_loss: 0.2815 - val_coeff_determination: 0.7069 - val_lr: 3.0518e-08
```

最後の方を出力しましたが、決定係数0.7くらいです。ハイパーパラメーターをいじればさらによいのかも。論文通りにSoftmax関数を最終層で使うとうまくいかないため、relu関数に変更するとうまくいきました。（なんでだろう？）

最後に

この技術での[分子生成技術](https://kivantium.net/deep-for-chem　"分子生成技術")も生まれたり、今後さらに応用が進むのではないかなと思っております。

参考

・SMILES2Vec: An Interpretable General-Purpose Deep Neural Network for Predicting Chemical Properties （元論文）
https://arxiv.org/abs/1712.02034
・今回参考にしたプログラム
https://github.com/Abdulk084/Smiles2vec/blob/master/smiles2vec.ipynb
・化合物でもDeep Learningがしたい！
https://kivantium.net/deep-for-chem