RDKitを使って化学構造式からlogSを予測する

2019.3.8追記
RDKitでQiitaを探していたら、ほぼ同じことをやっている記事を見つけた。
https://qiita.com/maskot1977/items/2e98d76a0782a6e89bde

やりたいこと

• 欲しい物性値になるような化学構造を知りたい
• 具体的には、官能基を変えたらどうなるかな～？を事前に知りたい
• 化学構造式を与えたら予測値を返してくれるモデルを作る

※RDKitのインストールはほかの記事を参考にして下さい。

問題設定

分子構造から溶解度を予測するモデルを作る

下準備

RDKitのサイトに溶解度(logS：100gの水にその物質が溶けた量[g]をwとしたとき、$\log_{10}w$ の値。無限に解けても上限は2)の情報のついたsdfファイル(solbility_test.sdf, solbility_train.sdf)があるので、ダウンロードして適当な場所に置いておく。以下ではpythonのファイルと同じ場所にあるとする。

モジュールのインポートなど。特にここでは説明はいらないだろう。

RDKit_solbility.py

from rdkit import Chem
from rdkit.Chem import Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors
from rdkit.Chem import PandasTools

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format='retina'
plt.style.use('seaborn-darkgrid')
c_cycle=("#ea5415","#005192","#ffdf00","#1d7a21","#88593b","#737061")

SDFファイルをDataFrameに取り込む

この記事の本体。

RDKit_solbility.py

def sdf_to_df(input_sdf):
    df_sdf = PandasTools.LoadSDF(input_sdf, molColName='Structure').set_index('NAME').drop(['ID','SOL_classification'], axis=1)
    names = df_sdf.index

    mols = [ mol for mol in Chem.SDMolSupplier(input_sdf)]
    descLists = [desc_name[0] for desc_name in Descriptors._descList]
    desc_calc = MoleculeDescriptors.MolecularDescriptorCalculator(descLists)
    data = [desc_calc.CalcDescriptors(mol) for mol in mols]
    df_desc = pd.DataFrame(data, columns=descLists, index=names)

    df = pd.concat([df_sdf, df_desc], axis=1)
    return df

1行目(df_sdf = )：LoadSDFでsdfファイルをDataFrameとして読み込む。読み込んだsdfファイルに複数の化学構造が含まれていればすべて読み込んでくれる。分子構造以外も全部読んでくれて、ここではlogSなども読み込んでいる。いらない情報はdrop()で落としている。必要に応じて。
4行目(mols = )：sdfファイルに含まれる分子のリスト
5行目(descLists = )：RDKitで出力できるdescriptorのリストを取得
6行目(desc_calc = )：descListsに含まれるdescriptorを計算してくれるインスタンス(でいいのかな？)
7行目(data = )：各分子ごとにdescriptorを計算する
8行目(df_desc)：descriptorの計算結果をDataFrameに変換
10行目：読み込んだsdfファイルのDataFrame(分子構造とlogSが含まれている)とdescriptorのDataFrameを結合

このDataFrameをdf.head()でチラ見すると、

こんな感じになる。

機械学習する

読み込んでしまえば、あとは通常の機械学習の手続きだけど参考までに書いておく。

SDFファイルの読み込み

RDKit_solvility.py

df = sdf_to_df('./solubility.train.sdf').dropna()
df_test = sdf_to_df('./solubility.test.sdf').dropna()

target_keys = ['SOL']

y_df = df.loc[:, target_keys].astype(float)
x_df = df.loc[:, [desc_name[0] for desc_name in Descriptors._descList]]

y_test = df_test.loc[:, target_keys].astype(float)
x_test = df_test.loc[:, [desc_name[0] for desc_name in Descriptors._descList]]

features = x_df.columns
samples  = x_df.index

x_df, y_dfにそれぞれ訓練データ、x_test, y_testはテストデータ。
手持ちのデータでわざわざ訓練データとテストデータを分けて保存しておくことはないと思うので、通常は読み込んだあとにtrain_test_splitなどを使って分けることになる。

回帰

RandomForestを使う。

RDKit_solvility.py

from sklearn.ensemble import RandomForestRegressor
reg = RandomForestRegressor(n_estimators=100, max_depth=500, random_state=0)
reg.fit(x_df, y_df.values[:,0])

回帰の本体。結果が芳しくなければパラメタや他の回帰方法を試せばいい。

可視化1

重要なDescriptorの可視化。全部のDescriptorを表示するようにするとえらいことになるので、
maxFeaturesかminImportanceで表示を制限している。

RDKit_solvility.py

from sklearn.metrics import r2_score

def plot_feature_importances(features, model, maxFeatures=-1, minImportance=-1, text=''):
    plt.figure(figsize=(16,10))
    plt.rcParams["font.size"] = 16
    
    df = pd.DataFrame({'features':features,
                       'importances':model.feature_importances_})
    df = df.set_index('features')
    #sort
    df = df.sort_values('importances', ascending=True)
    #normalize
    func_norm = lambda x: x/max(x)
    df = df.apply(func_norm)

    #minImportace
    df = df.applymap(lambda x:np.nan if abs(x)<minImportance else x).dropna(axis=0, how='all')
    
    if maxFeatures>0:
        df = df[-maxFeatures:]
    
    n_features = len(df.index)
    plt.barh(range(n_features), df['importances'], align='center')
    plt.yticks(np.arange(n_features), df.index)
    plt.xlabel(text + ' Feature importance (Random Forest)')
    plt.ylabel('Feature')

plot_feature_importances(features, reg, minImportance=0.01)

可視化2

observe-predictでプロット、$R^2$を表示する。

RDKit_solvility.py

ytests = []
ypreds = []
ytests = y_test.values[:,0]
ypreds = reg.predict(x_test)

plt.figure(figsize=(16,9))
plt.rcParams["font.size"] = 18

plt.xlabel('observe')
plt.ylabel('predict')
r2 = r2_score(np.array(ytests), np.array(ypreds))
plt.plot(ytests, ypreds, 'o', ms=15, color=c_cycle[0], alpha=0.3, label='$R^2$ = '+str(r2))
plt.legend()

結果

重要Descriptor

…**アカーーーン！！！** これだと「logSを予測するにはlogPが重要です」と言っているようなもの。 logSとlogPは負の相関があるのは**ほぼ**自明なので、構造から計算するlogPの精度が高いと**ほぼ同値**のことを言っていることになってしまう。まあ、「化学構造入れたらlogSを予測できるモデルを作る」というならそれでもいいんだけれども。

さすがに面白くないので、logPのデータを落としてもう一回。

RDKit_solvility.py

x_df = x_df.drop('MolLogP', axis=1)
x_test = x_test.drop('MolLogP', axis=1)

MR,分子屈折率？が重要？ 他には部分電荷や表面積の一種など。 何が来ても意外だけど、意外な結果であった。

R^2

おお、意外と高い。 でも、値として1違うと溶解量が1桁違うので、実用的なモデルにするには$R^2=0.999$くらい必要になるのではないだろうか…