これまでの記事では、Fingerprintの可視化や構造式描画のカスタマイズについてまとめてきた。
今回は、それらでまとめてきた手法を駆使して、機会学習による化合物の物性予測根拠の可視化を実施する。
具体的には、第43回ケモインフォマティクス討論会での講演「フィンガープリントによる化合物の予測根拠可視化（三井化学）」の内容を再現する。
講演（予稿）の内容
線形の回帰手法（機械学習を含む）を用いた化合物の物性予測においてその予測根拠を可視化する手法を提案している。具体的には"フィンガープリントを記述子として用いた機械学習モデルにおいて、各bitの寄与を分子構造の広がりを考慮して統合し化合物構造上に可視化する"というもの。
予稿中では、水溶解度のデータを検証用として、予測根拠の可視化手法としての有用性を化学的な知見で評価している。
講演予稿集より引用
解析条件など
データ
1290分子の水溶解度（LogS）データ（sdfファイル形式）
（金子研究室HPのデータセットを利用）
変数
radius=2 の Morganフィンガープリントを使用
前処理
- 以下の変数を削除
- 分散が0の変数
- 互いに相関係数が高い（0.95以上）変数ペアの内、片方の変数
（金子研モジュール使用） - 1090以上の分子で同じ値をとる変数
- 変数は全て標準化（平均：0, 分散：1）
- Borutaにより変数選択した（参考：Borutaによる変数選択）
解析手法
- テストデータの数：200/1290
- 予測器：線形サポートベクターマシン
- ハイパーパラメータ最適化：5-fold CV で実施
寄与の計算方法
ある原子 $i$ に対する部分構造の寄与を統合した値 $A_i$ を次の式で算出する。
$$ A_i = \sum_{n=1}^N(C_n\times\frac{1}{f_n}\times\frac{1}{x_n}) $$
$f$：部分構造を構成する原子に一つ以上 $i$ を含んでいるフィンガープリントのリスト
$C_n$：各フィンガープリントの寄与
$f_n$：分子中に含まれる各部分構造の数($n=1,2,...N$)
$x_n$：各部分構造に含まれる原子数
可視化方法
寄与の大きさに応じて以下の色で原子にハイライトする。
- 正に寄与している原子：赤
- 負に寄与している原子：青
前準備
- データの取得（金子研究室HPのデータセットを利用）
- モジュールのインストール
conda install -c conda-forge rdkit
pip install dcekit
pip install borutapyコード
データセットの準備〜前処理まで
import os
import pandas as pd
import numpy as np
import matplotlib.figure as figure
import matplotlib.pyplot as plt
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG
from boruta import BorutaPy
from sklearn import svm
from sklearn.preprocessing import StandardScaler
from sklearn import preprocessing
from sklearn.ensemble import RandomForestRegressor as rf
from sklearn.model_selection import train_test_split,GridSearchCV
from sklearn.metrics import r2_score, mean_absolute_error
from dcekit.variable_selection import search_high_rate_of_same_values, search_highly_correlated_variables
suppl = Chem.SDMolSupplier(os.getcwd() + '/../../data/logSdataset1290_2d.sdf')
mols_list = [mol for mol in suppl if mol is not None]
# ECFP4フィンガープリントの作成（Extended Connectivity Fingerprint）
ECFP4 = [AllChem.GetMorganFingerprintAsBitVect(mol, 2,1024) for mol in mols_list]
df = pd.DataFrame(np.array(ECFP4, int))
#溶解度(logS)の値を取得
df['logS'] = [mol.GetProp('logS') if 'logS' in mol.GetPropNames() else 'None' for mol in mols_list]
df = df.astype(float)
#前処理=====================
#欠損のある変数を削除
df.dropna(axis=1,inplace=True)
#分散が０の変数削除
del_num1 = np.where(df.var() == 0)
df = df.drop(df.columns[del_num1], axis=1)
#変数選択（互いに相関関係にある変数の内一方を削除）
threshold_of_r = 0.95 #変数選択するときの相関係数の絶対値の閾値
corr_var = search_highly_correlated_variables(df, threshold_of_r)
df.drop(df.columns[corr_var], axis=1, inplace=True)
#同じ値を多くもつ変数の削除
#1090化合物をテストデータとするため、1090化合物以上で同じ値をとる変数を削除
rate_of_same_value = []
threshold_of_rate_of_same_value = 1090
for col in df.columns:
same_value_number = df[col].value_counts()
rate_of_same_value.append(float(same_value_number[same_value_number.index[0]] / df.shape[0]))
del_var_num = np.where(np.array(rate_of_same_value) >= threshold_of_rate_of_same_value)
df.drop(df.columns[del_var_num], axis=1, inplace=True)
#前処理終わり=====================
#訓練データと学習データに分割
df_X = df.drop(['logS'],axis=1)
df_y = df.loc[:,['logS']]
X = df_X.values
y = df_y.values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)
#標準化
scaler_X = StandardScaler()
scaler_y = StandardScaler()
scaled_X_train = scaler_X.fit_transform(X_train)
scaled_y_train = scaler_y.fit_transform(y_train)
scaled_X_test = scaler_X.transform(X_test)
scaled_y_test = scaler_y.transform(y_test)borutaによる変数選択
（参考：Borutaによる変数選択）
#Borutaによる変数選択====================
#pパーセンタイルの最適化（金子研の方法）
corr_list = []
for n in range(10000):
shadow_features = np.random.rand(scaled_X_train.shape[0]).T
corr = np.corrcoef(scaled_X_train, shadow_features,rowvar=False)[-1]
corr = abs(corr[corr < 0.95])
corr_list.append(corr.max())
corr_array = np.array(corr_list)
perc = 100 * (1-corr_array.max())
#Borutaの実施
feat_selector = BorutaPy(rf(),
n_estimators='auto', # 特徴量の数に比例して、木の本数を増やす
verbose=0,
alpha=0.05, # 有意水準
max_iter=50, # 試行回数
perc = perc, #ランダム生成変数の重要度の何％を基準とするか
random_state=0
)
feat_selector.fit(scaled_X_train, scaled_y_train)
selected_bit = feat_selector.support_
selected_bit_num = df_X.columns[selected_bit]
scaled_X_train_br = scaled_X_train[:,selected_bit]
scaled_X_test_br = scaled_X_test[:,selected_bit]
print('pパーセンタイル:',round(perc,2))
print('boruta後の変数の数:',scaled_X_train_br.shape[1])
"""output
pパーセンタイル: 81.6
boruta後の変数の数: 136
"""予測モデル構築と精度検証
#クロスバリデーションによるハイパーパラメータ最適化
params = {'C':2 ** np.arange(-5, 5, 3, dtype=float),
'epsilon':2 ** np.arange(-10, 0, 3, dtype=float)
}
model_cv = GridSearchCV(svm.SVR(kernel='linear'), params,cv=5, n_jobs=-1)
model_cv.fit(scaled_X_train_br, scaled_y_train)
model = svm.SVR(kernel='linear',
C=model_cv.best_params_['C'],
epsilon=model_cv.best_params_['epsilon'])
#学習と予測精度の調査
model.fit(scaled_X_train_br, scaled_y_train)
scaled_pred_y_train = np.ndarray.flatten(model.predict(scaled_X_train_br))
scaled_pred_y_test = np.ndarray.flatten(model.predict(scaled_X_test_br))
pred_y_train = scaler_y.inverse_transform(scaled_pred_y_train)
pred_y_test = scaler_y.inverse_transform(scaled_pred_y_test)
print('r2_train:{}'.format(round(r2_score(y_train, pred_y_train),3)))
print('r2_test:{}'.format(round(r2_score(y_test, pred_y_test),3)))
print('MAE_train:{}'.format(round(mean_absolute_error(y_train, pred_y_train)3)))
print('MAE_test:{}'.format(round(mean_absolute_error(y_test, pred_y_test),3)))
plt.figure(figsize=(6,6))
plt.scatter(y_train, pred_y_train,color='blue',alpha=0.3,label='training data')
plt.scatter(y_test, pred_y_test, color='red',alpha=0.3,label='test data')
y_max_ = max(y_train.max(), pred_y_train.max(),y_test.max(), pred_y_test.max())
y_min_ = min(y_train.min(), pred_y_train.min(),y_test.min(), pred_y_test.min())
y_max = y_max_ + (y_max_ - y_min_) * 0.1
y_min = y_min_ - (y_max_ - y_min_) * 0.1
plt.plot([y_min , y_max],[y_min, y_max], 'k-')
plt.ylim(y_min, y_max)
plt.xlim(y_min, y_max)
plt.xlabel('logS(observed)',fontsize=20)
plt.ylabel('logS(predicted)',fontsize=20)
plt.legend(loc='best',fontsize=15)
plt.title('yyplot(logS)',fontsize=20)
plt.savefig('yyplot.png')
plt.show()予測精度は以下の通りで、予稿集と同程度になった。
| r2（決定係数） | MAE（誤差平均） | |
|---|---|---|
| train data | 0.729 | 0.84 |
| test data | 0.615 | 0.984 |
回帰係数の可視化
#全データを使ってモデルを再構築
scaler_X = StandardScaler()
scaler_y = StandardScaler()
scaled_X = scaler_X.fit_transform(X[:,selected_bit])
scaled_y = scaler_y.fit_transform(y)
model.fit(scaled_X, scaled_y)
pred_y = model.predict(scaled_X)
#回帰係数の可視化
selected_bit = feat_selector.support_
selected_bit_num = df_X.columns[selected_bit]
#回帰係数の取得
regression_coef = model.coef_[0]
#棒グラフの作成
x = np.arange(len(selected_bit_num))
width = 0.95
fig, ax = plt.subplots(figsize=(20,5))
bar = ax.bar(x, regression_coef, width)
ax.set_xticks(x)
ax.set_xticklabels(selected_bit_num,fontsize=5)
ax.set_ylabel('Regression Coefficient',fontsize=20)
ax.set_title('Regression Coefficient each bits',fontsize=25)
plt.grid(linestyle='--',axis='y')
plt.savefig('Regression Coefficient.png',bbox_inches='tight')
plt.show()
寄与の計算と可視化
寄与の計算（おさらい）
ある原子 $i$ に対する部分構造の寄与を統合した値 $A_i$ を次の式で算出する。
$$ A_i = \sum_{n=1}^N(C_n\times\frac{1}{f_n}\times\frac{1}{x_n}) $$
$f$：部分構造を構成する原子に一つ以上 $i$ を含んでいるフィンガープリントのリスト
$C_n$：各フィンガープリントの寄与
$f_n$：分子中に含まれる各部分構造の数($n=1,2,...N$)
$x_n$：各部分構造に含まれる原子数
可視化方法（おさらい）
寄与の大きさに応じて以下の色で原子にハイライトする。
- 正に寄与している原子：赤
- 負に寄与している原子：青
化合物の描画のカスタマイズについては、以下記事を参照
（rdMolDraw2Dモジュールを使って構造式描画をカスタマイズ）
可視化
適当に一分子ピックアップして可視化してみる。
#ビットと回帰係数の辞書を作っておく
bit_coef = dict(zip(selected_bit_num, regression_coef))
#適当な分子をピックアップ
mol = mols_list[1000]
print(Chem.MolToSmiles(mol))
#フィンガープリントを計算
bitI_morgan = {}
fp_morgan = AllChem.GetMorganFingerprintAsBitVect(mol, 2, 1024, bitInfo=bitI_morgan)
#寄与のあるビットの抜き出し
bit_list = list(set(selected_bit_num)&set(bitI_morgan.keys()))
print('寄与のあるビットの数：',len(bit_list))
#寄与を格納する配列を生成
Ai_list = np.zeros(mol.GetNumAtoms())
for bit in bit_list:
#フィンガープリントの寄与
Cn = bit_coef[bit]
#分子中に含まれる部分構造の数
fn = len(bitI_morgan[bit])
for part in bitI_morgan[bit]:
#fingerprintのradiusが0の時はfingerprintの中心原子を抜き出し、寄与を加算する。
if part[1]==0:
i = part[0]
xn = 1
Ai_list[i] += Cn / fn / xn
#fingerprintのradiouが0以上の時は、該当する原子のリストを抜き出し、寄与を加算する
else:
amap={}
env = Chem.FindAtomEnvironmentOfRadiusN(mol,
radius=part[1],
rootedAtAtom=part[0])
submol=Chem.PathToSubmol(mol,env,atomMap=amap)
#各部分構造に含まれる原子数
xn = len(list(amap.keys()))
for i in amap.keys():
Ai_list[i] += Cn / fn / xn
#正規化っぽいことをする
#絶対値の最大値が0.5になるように調整する
Ai_list = Ai_list / abs(Ai_list).max() * 0.5
print('寄与の値')
print(Ai_list)
# ハイライトする原子と色，円の半径を設定
atoms = [i for i in range(len(Ai_list))]
atom_colors = dict()
for i in atoms:
if Ai_list[i] > 0:
atom_colors[i] = (1,1-Ai_list[i],1-Ai_list[i])
else:
atom_colors[i] = (1+Ai_list[i],1+Ai_list[i],1)
# コンテナの準備
view = rdMolDraw2D.MolDraw2DSVG(300,350)
tm = rdMolDraw2D.PrepareMolForDrawing(mol)
view.DrawMolecule(tm,
highlightAtoms=atoms,
highlightAtomColors=atom_colors,
highlightBonds=[],
highlightBondColors={})
#ファイナライズ、保存、描画
view.FinishDrawing()
svg = view.GetDrawingText()
with open('highlighted_sample.svg', 'w') as f:
f.write(svg)
SVG(svg)結果
化学的な知見とも合うので可視化手法として有用と言える。
寄与の正負を問わないのであればランダムフォレスト系の手法にも適用できそう。
参考
- Borutaによる変数選択
- rdMolDraw2Dモジュールを使って構造式描画をカスタマイズ
- フィンガープリントによる化合物の予測根拠可視化
- RDKitでの構造式描画を詳しく解説
- RDkitで分子フィンガープリントの生成
関連書籍
リンク
リンク
リンク