以前、構造物性相関解析（QSPR）でフィンガープリントを変数として回帰分析した際に、化学構造のどの部分が物性に寄与しているかを可視化するという記事を書いた（フィンガープリントによる化合物の予測根拠可視化）。
ただ、以前用いた方法は線形回帰における回帰係数をもとに、各部分構造が正/負に寄与しているかを赤/青で表すものなので、回帰係数の算出できない非線形手法には使えなかった。
そこで、SHAP値を使うことで、非線形の手法であっても化学構造の寄与が可視化できるか検討した。
SHAPを使うと、ある機械学習モデルが導き出した予測値に対して、それぞれの特徴量の影響をどれだけ受けたか（寄与の度合い）を求めることができる。
（詳しくはSHAP を用いた機械学習への解釈性付与）
解析条件など
データ
1290分子の水溶解度（LogS）データ（sdfファイル形式）
（金子研究室HPのデータセットを利用）
変数
radius=2 の Morganフィンガープリントを使用
前処理
- 以下の変数を削除
- 分散が0の変数
- 互いに相関係数が高い（0.95以上）変数ペアの内、片方の変数（金子研モジュール使用）
解析手法
- 精度検証方法：Hold-out法（テストデータ比率：0.25）
- 予測器：XGBoost
- ハイパーパラメータ最適化：なし
寄与の計算方法
ある原子 $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$：各フィンガープリントの寄与（SHAPによって算出した値）
- $f_n$：分子中に含まれる各部分構造の数($n=1,2,...N$)
- $x_n$：各部分構造に含まれる原子数
可視化方法
寄与の大きさに応じて以下の色で原子にハイライトする。
- 正に寄与している原子：赤
- 負に寄与している原子：青
モジュールインポート〜精度検証まで
import shap
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
import xgboost
from sklearn import svm
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_absolute_error
from dcekit.variable_selection import search_highly_correlated_variables
# データの読み込み
df = pd.read_csv(r'./logSdataset1290.csv',index_col=0)
# フィンガープリント生成
mols_list = [Chem.MolFromSmiles(smiles) for smiles in df.index]
ECFP4 = [AllChem.GetMorganFingerprintAsBitVect(mol, 2,1024) for mol in mols_list]
df_ECFP4 = pd.DataFrame(np.array(ECFP4, int), index=df.index)
df = pd.concat([df_ECFP4,df['logS']],axis=1)
#欠損のある変数、分散が0の変数を削除
df.dropna(axis=1,inplace=True)
del_num1 = np.where(df.var() == 0)
df = df.drop(df.columns[del_num1], axis=1)
#変数選択（互いに相関係数が0.95以上の変数の内一方を削除）
corr_var = search_highly_correlated_variables(df, 0.95)
df.drop(df.columns[corr_var], 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.25, random_state=0)
#学習と予測精度の調査
model = xgboost.XGBRegressor()
model.fit(X_train, y_train)
pred_y_train = np.ndarray.flatten(model.predict(X_train))
pred_y_test = np.ndarray.flatten(model.predict(X_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_train:0.969
r2_test:0.717
MAE_train:0.232
MAE_test:0.841
十分な精度だと判断して次に進む。
SHAP値を算出する
全てのデータを使ってモデルを再構築し、SHAP値を計算する。
#全データを使ってモデルを再構築
model = xgboost.XGBRegressor()
model.fit(df_X, df_y)
# TreeExplainer でモデルと学習に使ったデータを渡してオブジェクトを作る。
explainer = shap.TreeExplainer(model=model,
feature_perturbation='interventional',
model_output='margin')
# shap_valueを計算する。
shap_values = explainer.shap_values(df_X)force_plot を描いてみる。
force_plotでは、何も特徴量の値が分からない状態の予測値（base value）から、各特徴量の寄与によって予測値がどう動いたかが視覚的に分かる
この下記のサンプルに関しては、ベースとなる溶解度の予測値（-2.728）に対して、ビット番号90や726などの寄与によって正に、ビット番号356や265などの寄与によって腑に負に予測値が押し上げ（押し下げ）られて最終的に-1.43という予測値となったことが読み取れる。
#index : 何番目のサンプルの情報をplotするか指定
index = 1000
shap.force_plot(base_value=explainer.expected_value,
shap_values=shap_values[index],
features=df_X.iloc[index,:],
)
SHAP値を化学構造に反映させる
#index : 何番目のサンプルの情報をplotするか指定
index = 1000
mol = Chem.MolFromSmiles(df.index[index])
print(df.index[index])
# ビットとSHAP値の辞書を作っておく
selected_bit_num = df_X.columns
bit_coef = dict(zip(selected_bit_num, shap_values[index]))
#フィンガープリントを計算
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)
Cc1c[nH]c(=O)[nH]c1=O
寄与のあるビットの数： 19
寄与の値
[-0.06748477 -0.33492273 -0.17817281 0.5 -0.32354564 0.32335094
0.31972844 -0.43044103 0.19928591]
極性の高い部分構造は、水への溶解度を高める方向、極性の低い炭素の部分は極性を下げる方向に寄与していることが読み取れるため、結果は妥当そう。
他にも色々可視化してみた。
以上
参考
参考HP
- SHAP を用いた機械学習への解釈性付与
- 技術レポート再現（フィンガープリントによる化合物の予測根拠可視化 | 三井化学）
- rdMolDraw2Dモジュールを使って構造式描画をカスタマイズ
- フィンガープリントによる化合物の予測根拠可視化
関連書籍
リンク
リンク
リンク
