Chemistry: Free Chemception models with RDKit and Keras

Garrett Goh及其同事發表兩篇有趣的論文https://arxiv.org/abs/1706.06689https://arxiv.org/abs/1706.06689和https://arxiv.org/abs/1710.02238，其中作者使用編碼圖像將分子讀入這些深度神經網絡 專為圖像識別而設計。 他們的代碼無法獲得，但用RDKit和Keras重建他們的Chemception網絡應該不會太難。 以下內容可能有點技術性和冗長，但最後有一些漂亮的圖像結果
首先是一些Python函式庫的導入。

1 2 3 4 5 6 7 8 9

from __future__ import print_function import rdkit from rdkit import Chem from rdkit.Chem import AllChem import pandas as pd import numpy as np from matplotlib import pyplot as plt %matplotlib inline print("RDKit: %s"%rdkit.__version__)

接下來是一些Keras的模組

1 2 3 4 5 6 7 8

import keras 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.preprocessing.image import ImageDataGenerator from keras.callbacks import ReduceLROnPlateau print("Keras: %s"%keras.__version__)

Preparing the Dataset

它在名為“smiles”的列中具有SMILE字串，並且在名為PC_uM_values的列中具有DHFR抑製劑的IC50數據集的分子數據集測定值。 在讀取數據集後，使用pandas dataframes apply函數將RDKit分子添加到名為“mol”的列中。

1 2	data = pd.read_csv("Sutherland.csv") data["mol"] = data["smiles"].apply(Chem.MolFromSmiles)

ps. Sutherland_DHFR.csv檔案

Chemception以Inception模塊命名，它將用於神經網絡。第一步是將分子編碼成“圖像”。下面的函數採用RDKit mol並將分子圖編碼為具有4個通道的圖像。

在讀取分子後，計算Gasteiger電荷併計算2D繪圖坐標。它們通常在生成分子的描述之前進行計算，但我們需要它們“原始”，因此它們被提取到coords矩陣。

vect矩陣被定義為三維陣列(image_width，image_height，4)並設定為零矩陣。然後每層編碼不同的分子資訊。channel的第0層填入有關於bonds的資訊並編碼bondorder。接下來的三層用原子序，Gasteiger電荷和hybridization來編碼。選擇這些資訊的原因是這種組合似乎在大多數測試中都相當順利，但當然也可以嘗試其他組合。

當我們有一個4通道圖像時，可以使用Keras的標準圖像處理工具來獲得適當的訓練。當然可以添加更多頻道，但這需要對Keras模塊進行一些額外的處理。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

def chemcepterize_mol(mol, embed=20.0, res=0.5):     dims = int(embed*2/res)     cmol = Chem.Mol(mol.ToBinary())     cmol.ComputeGasteigerCharges()     AllChem.Compute2DCoords(cmol)     coords = cmol.GetConformer(0).GetPositions()     vect = np.zeros((dims,dims,4))     #Bonds first     for i,bond in enumerate(mol.GetBonds()):         bondorder = bond.GetBondTypeAsDouble()         bidx = bond.GetBeginAtomIdx()         eidx = bond.GetEndAtomIdx()         bcoords = coords[bidx]         ecoords = coords[eidx]         frac = np.linspace(0,1,int(1/res*2)) #         for f in frac:             c = (f*bcoords + (1-f)*ecoords)             idx = int(round((c[0] + embed)/res))             idy = int(round((c[1]+ embed)/res))             #Save in the vector first channel             vect[ idx , idy ,0] = bondorder     #Atom Layers     for i,atom in enumerate(cmol.GetAtoms()):             idx = int(round((coords[i][0] + embed)/res))             idy = int(round((coords[i][1]+ embed)/res))             #Atomic number             vect[ idx , idy, 1] = atom.GetAtomicNum()             #Gasteiger Charges             charge = atom.GetProp("_GasteigerCharge")             vect[ idx , idy, 3] = charge             #Hybridization             hyptype = atom.GetHybridization().real             vect[ idx , idy, 2] = hyptype     return vect

為了更好地理解代碼的作用，我們嘗試“ chemcepterize”一個分子並將其顯示為圖像。 embedding和resolution設定得低於最終數據集。 Matplotlib僅支援RGB，因此僅使用前三個通道。

1 2 3 4

mol = data["mol"][0] v = chemcepterize_mol(mol, embed=10, res=0.2) print(v.shape) plt.imshow(v[:,:,:3])

我們可以看到不同的原子和鍵的顏色不同，因為它們具有不同的化學資訊。 下一步是“chemcepterize”整個RDKit分子集合，並在dataframe中添加帶有“ images”的新column

1 2 3

def vectorize(mol):     return chemcepterize_mol(mol, embed=12) data["molimage"] = data["mol"].apply(vectorize)

將數據集分化成訓練跟測試集而且最後numpy陣列的形狀（樣本，高度，寬度，通道）

seq = np.ones(756)
for i in range(147):
    seq[i] = 0
np.random.shuffle(seq)
for i in range(756):
    data.at[i,'split'] = seq[i]
(pandas的column增加split來label 0代表testing data跟1代表training data)

1 2 3 4

X_train = np.array(list(data["molimage"][data["split"]==1])) X_test = np.array(list(data["molimage"][data["split"]==0])) print(X_train.shape) print(X_test.shape)

我們還需要準備要預測DHFR抑製劑的IC50值。 數據被轉換為對數空間(log space)，來自scikit-learn的強大的縮放器用於將數據縮放到-1到1之間(像這樣範圍的神經網絡，訓練會更容易)。

1 2 3 4 5 6 7 8 9

assay = "PC_uM_value" y_train = data[assay][data["split"]==1].values.reshape(-1,1) y_test = data[assay][data["split"]==0].values.reshape(-1,1) from sklearn.preprocessing import RobustScaler rbs = RobustScaler(with_centering=True, with_scaling=True, quantile_range=(5.0, 95.0), copy=True) y_train_s = rbs.fit_transform(np.log(y_train)) y_test_s = rbs.transform(np.log(y_test)) h = plt.hist(y_train_s, bins=20) plt.show()

Building the Chemception model
現在我們將建立chemception的有趣部分。為此，我們將使用Inception模塊。這些都是有趣的架構。模塊不是選擇/優化卷積內核大小，而是將數據流分成三個。第一個塔是圖像通道/先前內核與捲積1的重組，而另外兩個分別具有3,3和5,5的內核大小。第一個塔與MaxPooling2D層組合，而其他塔與先前的1個1卷積層組合。藉由重新組合先前層的特徵到較少數量的通道來縮小dimension。

第一個初始模塊是特殊的，因為我已經刪除了MaxPooling層。理由是編碼的特徵不是自然排序的。正電荷是否比負電荷更有意義？氧氣是否比碳和氮更重要，因為它具有更高的原子量？我想不是，所以在開始使用帶有max pooling的標准單元之前，利用第一層重新組合通道的方式讓網絡來理解。

在一些輔助函數中定義Inception模塊之後，使用Keras功能API構建網絡，其中定義了Keras張量對象(tensor objects)的流。

在初始模塊之後，我在初始模塊中放置了一個跨越每個內核輸出的max pooling layer，希望在不同的內核中創建一些特徵。我們感興趣的是特徵是否存在，但可能不是圖像的確切位置。將輸出平坦化為一維，並將其輸入標準密集網絡，其中具有RELU激活函數的100個神經元，和具有線性激活函數的單一個輸出神經元，因為我們的目標是回歸模型。

最後，計算圖與模型類一起使用，其中定義了輸入和輸出。

1 2	input_shape = X_train.shape[1:] print(input_shape)

1 2 3 4 5 6 7 8

def Inception0(input):     tower_1 = Conv2D(16, (1, 1), padding='same', activation='relu')(input)     tower_1 = Conv2D(16, (3, 3), padding='same', activation='relu')(tower_1)     tower_2 = Conv2D(16, (1, 1), padding='same', activation='relu')(input)     tower_2 = Conv2D(16, (5, 5), padding='same', activation='relu')(tower_2)     tower_3 = Conv2D(16, (1, 1), padding='same', activation='relu')(input)     output = keras.layers.concatenate([tower_1, tower_2, tower_3], axis=-1)     return output

1 2 3 4 5 6 7 8 9

def Inception(input):     tower_1 = Conv2D(16, (1, 1), padding='same', activation='relu')(input)     tower_1 = Conv2D(16, (3, 3), padding='same', activation='relu')(tower_1)     tower_2 = Conv2D(16, (1, 1), padding='same', activation='relu')(input)     tower_2 = Conv2D(16, (5, 5), padding='same', activation='relu')(tower_2)     tower_3 = MaxPooling2D((3, 3), strides=(1, 1), padding='same')(input)     tower_3 = Conv2D(16, (1, 1), padding='same', activation='relu')(tower_3)     output = keras.layers.concatenate([tower_1, tower_2, tower_3], axis=-1)     return output

1 2 3 4 5 6 7 8 9 10 11

input_img = Input(shape=input_shape) x = Inception0(input_img) x = Inception(x) x = Inception(x) od=int(x.shape[1]) x = MaxPooling2D(pool_size=(od,od), strides=(1,1))(x) x = Flatten()(x) x = Dense(100, activation='relu')(x) output = Dense(1, activation='linear')(x) model = Model(inputs=input_img, outputs=output) print(model.summary())

__________________________________________________________________________________________________
 Layer (type) Output Shape Param # Connected to
 ==================================================================================================
 input_1 (InputLayer) (None, 48, 48, 4) 0
 __________________________________________________________________________________________________
 conv2d_1 (Conv2D) (None, 48, 48, 16) 80 input_1[0][0]
 __________________________________________________________________________________________________
 conv2d_3 (Conv2D) (None, 48, 48, 16) 80 input_1[0][0]
 __________________________________________________________________________________________________
 conv2d_2 (Conv2D) (None, 48, 48, 16) 2320 conv2d_1[0][0]
 __________________________________________________________________________________________________
 conv2d_4 (Conv2D) (None, 48, 48, 16) 6416 conv2d_3[0][0]
 __________________________________________________________________________________________________
 conv2d_5 (Conv2D) (None, 48, 48, 16) 80 input_1[0][0]
 __________________________________________________________________________________________________
 concatenate_1 (Concatenate) (None, 48, 48, 48) 0 conv2d_2[0][0]
 conv2d_4[0][0]
 conv2d_5[0][0]
 __________________________________________________________________________________________________
 conv2d_6 (Conv2D) (None, 48, 48, 16) 784 concatenate_1[0][0]
 __________________________________________________________________________________________________
 conv2d_8 (Conv2D) (None, 48, 48, 16) 784 concatenate_1[0][0]
 __________________________________________________________________________________________________
 max_pooling2d_1 (MaxPooling2D) (None, 48, 48, 48) 0 concatenate_1[0][0]
 __________________________________________________________________________________________________
 conv2d_7 (Conv2D) (None, 48, 48, 16) 2320 conv2d_6[0][0]
 __________________________________________________________________________________________________
 conv2d_9 (Conv2D) (None, 48, 48, 16) 6416 conv2d_8[0][0]
 __________________________________________________________________________________________________
 conv2d_10 (Conv2D) (None, 48, 48, 16) 784 max_pooling2d_1[0][0]
 __________________________________________________________________________________________________
 concatenate_2 (Concatenate) (None, 48, 48, 48) 0 conv2d_7[0][0]
 conv2d_9[0][0]
 conv2d_10[0][0]
 __________________________________________________________________________________________________
 conv2d_11 (Conv2D) (None, 48, 48, 16) 784 concatenate_2[0][0]
 __________________________________________________________________________________________________
 conv2d_13 (Conv2D) (None, 48, 48, 16) 784 concatenate_2[0][0]
 __________________________________________________________________________________________________
 max_pooling2d_2 (MaxPooling2D) (None, 48, 48, 48) 0 concatenate_2[0][0]
 __________________________________________________________________________________________________
 conv2d_12 (Conv2D) (None, 48, 48, 16) 2320 conv2d_11[0][0]
 __________________________________________________________________________________________________
 conv2d_14 (Conv2D) (None, 48, 48, 16) 6416 conv2d_13[0][0]
 __________________________________________________________________________________________________
 conv2d_15 (Conv2D) (None, 48, 48, 16) 784 max_pooling2d_2[0][0]
 __________________________________________________________________________________________________
 concatenate_3 (Concatenate) (None, 48, 48, 48) 0 conv2d_12[0][0]
 conv2d_14[0][0]
 conv2d_15[0][0]
 __________________________________________________________________________________________________
 max_pooling2d_3 (MaxPooling2D) (None, 1, 1, 48) 0 concatenate_3[0][0]
 __________________________________________________________________________________________________
 flatten_1 (Flatten) (None, 48) 0 max_pooling2d_3[0][0]
 __________________________________________________________________________________________________
 dense_1 (Dense) (None, 100) 4900 flatten_1[0][0]
 __________________________________________________________________________________________________
 dense_2 (Dense) (None, 1) 101 dense_1[0][0]
 ==================================================================================================
 Total params: 36,153
 Trainable params: 36,153
 Non-trainable params: 0
 __________________________________________________________________________________________________ 

最後，模型最終沒有很多參數，只有36.153，比我在更小的數據集上訓練的其他參數要低很多。 這是因為考慮到數據集相當小的關係，所以將每個塔的內核數量從64減少到16。 這些圖像可能比貓的照片簡單得多，所以無論如何它都可以工作。

為了優化，使用Adam優化器和均方誤差作為損失函數。

1 2	optimizer = Adam(lr=0.00025) model.compile(loss="mse", optimizer=optimizer)

下一部分對於避免過度擬合至關重要。這裡，ImageDataGenerator物件用於在訓練之前執行圖像的隨機旋轉和翻轉，作為增強訓練數據集的一種方式。藉由這個方式，網絡將學習如何處理旋轉並且看到不同方向的特徵將有助於模型更好地歸納。不使用增強訓練數據集將導致完全overfit模型。我們沒有在圖像中編碼立體化學訊息，否則翻轉應該藉由其他方式達成。

訓練集連接到50倍的長度，以具有合理的epochs大小。

1 2 3 4 5 6 7 8 9 10 11 12

from image import ImageDataGenerator generator = ImageDataGenerator(rotation_range=180,                                width_shift_range=0.1,height_shift_range=0.1,                                fill_mode="constant",cval = 0,                                horizontal_flip=True, vertical_flip=True,data_format='channels_last',                                ) #Concatenate for longer epochs Xt = np.concatenate([X_train]*50, axis=0) yt = np.concatenate([y_train_s]*50, axis=0) batch_size=128 g = generator.flow(Xt, yt, batch_size=batch_size, shuffle=True) steps_per_epoch = 10000/batch_size

Training the Chemception model
現在有趣的部分：訓練。 為了降低學習率，一旦驗證丟失開始穩定，我就可以使用ReduceLROnPlateau回調。

1 2 3 4 5 6

reduce_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.5,patience=10, min_lr=1e-6, verbose=1) history = model.fit_generator(g,                               steps_per_epoch=len(Xt)//batch_size,                               epochs=150,                               validation_data=(X_test,y_test_s),                               callbacks=[reduce_lr])

Epoch 1/150
 237/237 [==============================] - 31s 131ms/step - loss: 0.1241 - val_loss: 0.1208
 Epoch 2/150
 237/237 [==============================] - 28s 119ms/step - loss: 0.1135 - val_loss: 0.1068
 Epoch 3/150
 237/237 [==============================] - 28s 119ms/step - loss: 0.1123 - val_loss: 0.1071
 Epoch 4/150
.....
237/237 [==============================] - 29s 122ms/step - loss: 0.0362 - val_loss: 0.0361
 Epoch 150/150
 237/237 [==============================] - 29s 122ms/step - loss: 0.0355 - val_loss: 0.0360

可以保存和加載模型。 歷史物件被保存。

1 2 3 4 5 6 7

name = "Chemception_std_notebook_demo" model.save("%s.h5"%name) hist = history.history import pickle pickle.dump(hist, file("%s_history.pickle"%name,"w")) #from keras.model import load_model #model = load_model("%s.h5"%name)

可以從學習過程的圖表來判斷訓練的收斂性。 有些不尋常，當沒有正規化時：驗證損失在損失之前下降。 驗證集未被增強，因此由一些“完美”圖片組成，而可能需要網絡更長時間來處理所有旋轉而引入一些像素偽像。

1 2 3 4 5 6

for label in ['val_loss','loss']:     plt.plot(hist[label], label = label) plt.legend() plt.yscale("log") plt.xlabel("Epochs") plt.ylabel("Loss/lr")

Plotting and Evaluating the Performance

1 2 3 4 5 6 7 8

y_pred_t = rbs.inverse_transform(model.predict(X_train)) y_pred = rbs.inverse_transform(model.predict(X_test)) plt.scatter(np.log(y_train), y_pred_t, label="Train") plt.scatter(np.log(y_test), y_pred, label="Test") plt.xlabel("log(PC_uM)") plt.ylabel("predicted") plt.plot([-10,6],[-10,6]) plt.legend()

1 2 3 4

corr2 = np.corrcoef(np.log(y_test).reshape(1,-1), y_pred.reshape(1,-1))[0][1]**2 rmse = np.mean((np.log(y_test) - y_pred)**2)**0.5 print("R2 : %0.2F"%corr2) print("RMSE : %0.2F"%rmse)

R2 : 0.67
RMSE : 1.50

這看起來非常類似於之前使用其他方法獲得的此數據集的性能。

Visualizing the Layers
嘗試理解模型如何“看到”分子可能會很有趣。 為此，將舉例一個分子並繪製來自不同層的一些輸出。 在數據集中採用了最低IC50，數字143的化合物。

1 2 3

molnum = 143 molimage = np.array(list(data["molimage"][molnum:molnum+1])) mol = data["mol"][molnum]

分子看起來像這樣

1 2	from rdkit.Chem import Draw Draw.MolToImage(mol)

並且這個“chemcepterized”圖像如下所示

1	plt.imshow(molimage[0,:,:,:3])

第一個例子是第三層，它是1,1卷積餵給塔2中的3,3卷積層。

1 2 3 4 5 6 7 8 9

layer1_model = Model(inputs=model.input,                     outputs=model.layers[2].output) kernels1 = layer1_model.predict(molimage)[0] def plot_kernels(kernels):     fig, axes = plt.subplots(2,3, figsize=(12,8))     for i,ax in enumerate(axes.flatten()):         ax.matshow(kernels[:,:,i])         ax.set_title("Kernel %s"%i) plot_kernels(kernels1)

不同的內核（這裡是16個中的前6個）已經在Chemception圖像中重新組合了信息。 內核5似乎專注於bond，因為它在其他層中存在原子時已經去除了bond資訊。 內核1關注原子，對脂肪族碳大量被激活。 內核4最多與氯原子一起離開，但也含有bond資訊。 內核2和3似乎是空的。 也許它們是被其他分子激活，在當前分子中特徵不存在，或者它們可能是不需要的。 讓我們更深入......

1 2 3 4

for layer in [7,13,15,19,20]:     print("Layer %i"%layer)     plot_kernels(Model(inputs=model.input,outputs=model.layers[layer].output).predict(molimage)[0])     plt.show()

Layer 7


Layer 13


Layer 15


Layer 19


Layer 20


透過各層越來越抽象和具體的特徵產生一個大體上的趨勢。 氯原子似乎在許多內核中點亮，但其他人則關注其他特徵。 第19層中的內核0到2似乎專注於所有不是氯原子。 同一層中的內核5在雙鍵氧原子附近激活。 在總核心max pooling之前的最後一層似乎只關注在非常特定部分的分子。 核0可以是酰胺氧。 內核2和5的氯原子。 內核4似乎喜歡雙鍵氧，但僅來自羧酸基團，而不是酰胺。 內核3在末端羧酸的OH附近被激活。 這些只是從圖像中提取前6個特徵並前饋到密集層。

Summary
在這篇文章中，展示如何用幾行python代碼創建一個化學模型，使用RDKit進行化學反應，使用Keras進行神經網絡。 Sutherland DHFR數據集非常依賴分子類別，並且隨機劃分類別到訓練和測試集中並不顯示這些類型的化學模型是否可以承載SAR從一個類到另一個類的任何轉移，或者他們只是識別複合類並指定平均IC50。



參考
https://www.wildcardconsulting.dk/learn-how-to-teach-your-computer-to-see-chemistry-free-chemception-models-with-rdkit-and-keras/