Smiles2vecで物性予測をしよう。

はじめに

今回はSmiles2vecについて簡単に書かせていただきます。お手柔らかに。。。アドベントカレンダーにどこか登録したかった、、、
ご意見ご感想いただければ幸いです。

Smiles2vecとは？

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

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

化学構造の情報を持ったベクトルを入力変数として、物性を目的変数として予測という感じでしょうか？
それでは、やっていきます。

Smiles2vecの構造について

上でも話したようにNLPの分野の技術で、Seq2Seqという名前の技術です。無学ながら簡単に説明すると「機械対話や機械翻訳などのモデル」でよく使用されているそうです。この技術は再帰的ニューラルネットワークの考えに基づいたLSTMやGRUなどの層を使います。下図は元論文のSmiles2vecの構造です。この論文では、化合物セットと物性によって層の中身を変えて試行しているようです。

実際にやってみた。

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

・importとデータのセット

smiles2vec.py

from __future__ import print_function
import keras
from keras import backend as K
from keras.callbacks import ModelCheckpoint
from keras import optimizers
from keras.layers import Dense
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
from sklearn.metrics import roc_auc_score, roc_curve
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
import numpy as np
import pandas as pd
import scipy
import tensorflow as tf
import xlsxwriter
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 sklearn.utils import shuffle
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

smiles	Activity	split
OC(=O)C1=C(Cl)C=CC=C1Cl	2	1
OC(=O)C1=C(Cl)C=C(Cl)C=C1Cl	3	1
OC(=O)C1=CC=CC(=C1Cl)Cl	3	1
OC(=O)C1=CC(=CC=C1Cl)Cl	3	1
OC1=C(C=C(C=C1)N+=O)N+=O	4	1
OC(=O)C1=CC(=CC(=C1)Cl)Cl	3	1
OC(=O)C1=CC=C(Cl)C(=C1)Cl	3	1
NC(=O)C1=CC=CC=C1	2	1
ClC(Cl)(Cl)Cl	3	1
C=CCN=C=S	5	1
CC(O)CO	0	1
O=C1CCO1	3	1
CC1=CC(=O)C2=C(C=CC=C2)C1=O	5	1
CC1=C(Cl)C=CC(=C1)O	4	1
OCCC1=CC=CC=C1	2	1
CC(N)=O	1	1
NC1=CC=CC=C1	3	1
NC(N)=S	2	1
CCO	1	1
CCOS(=O)OCC	2	1
NC(=O)C1=CC=CC=C1O	3	1

結果

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関数に変更するとうまくいきました。（なんでだろう？）

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

最後に

いかかだったでしょうか？正直NLP分野は初めてだったのですが、面白かったので、参考書を買って勉強中です。この技術での分子生成技術も生まれたり、今後さらに応用が進むのではないかなと思っております。
また、今回の記事が面白ければ、下記の記事はわかりやすく、興味深い内容が書かれているので是非。
https://qiita.com/shionhonda/items/a748840b1a32e4812129
（この記事のおかげでいろいろと捗りました。この場を借りて感謝いたします。）

追記
いろいろ参考書をみたのですが、最終層がSoftmax関数ではないのは、encodeはするが、decodeする必要がないからですね！

参考

・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