mordred記述子を活用したケモインフォマティクス入門

今回の記事は下記のような読者に向けて書いています。

「pythonを活用してケモインフォマティクスを実践したい」
「mordred記述子の使い方を知りたい」
「mordred記述子の欠損値処理を詳しく知りたい」
「mordred記述子を使ったケモインフォマティクスの実践例を見てみたい」

個人的にmordred記述子を愛用してまして、化学構造をmordred記述子でベクトル化し、LightGBMでモデル構築というやり方が最も予測精度が出るのではないかと思っています。

一方で、mordred記述子は使い方が少し難しいので、今回解説記事を書くことにしました。

この記事でmordred記述子の使い方を理解して、是非使ってみてください。
早速解説します。

mordred記述子とは

mordred記述子は大阪大学が開発した分子記述子の一つです。
特徴としては、化学構造を高次元のベクトルに変換できます。

rdkit記述子 208次元
mordred記述子　1826次元

mordred記述子は化学構造を高次元のベクトルに変換できるので、より詳細に化学構造の特徴を表現できていると考えていいです。

次に使い方について解説していきます。

mordred記述子の使い方

使い方は、rdkit記述子と同様に「MOLファイル→ベクトル」に変換する際に使用します。

mordred記述子はanacondaにデフォルトで入っていないので、別途インストールが必要です。

mordred記述子のインストール

mordred記述子はrdkitと一緒に使います。（SMILES→MOLファイルへの変換はrdkitを活用）
なので下記コードを使って、rdkitがインストールされているpython環境にインストールしてください。

conda install conda-forge::mordred

またインストールできない場合は下記サイトを参考にしてください。

Mordred :: Anaconda.org

mordred記述子をインストールできたと思うので、実際にmordred記述子を使ってケモインフォマティクスを実践していきます。

mordred記述子を使ったケモインフォマティクスの実践

今回は、mordred記述子を使って水溶解度予測モデルを構築します。

目的変数を水溶解度、説明変数を化学構造（ベクトル）として、その２つを結ぶFを機械学習で構築します。

水溶解度 = F（化学構造）

目的変数の水溶解度は、水溶解度をS（mol/L)とすると、対数変換したlog₁₀Sの形で使用します。

データの読み込みからMOLファイルに変換まで

データの読み込みから、SMILES→MOLファイルへの変換は下記記事で解説してるので、是非ご覧ください。

rdkitを活用したケモインフォマティクス入門

今回の記事は下記のような読者に向けて書いています。「pythonを活用してケモインフォマティクスを実践したい」「rdkitで何ができるか知りたい」「rdkitを使ったケモインフォマティクスの実施例を見てみたい」 pythonを活用した...

念のため、下記にMOLファイルに変換するところまでpythonコードを記載しておきます。
（内容理解している方はコピペして使ってください）

import urllib.request # URLで指定したファイルをダウンロードするライブラリ
import pandas as pd
import numpy as np
from rdkit import Chem

url = 'http://modem.ucsd.edu/adme/data/databases/logS/data_set.dat'
urllib.request.urlretrieve(url, 'water_solubility.txt')
df = pd.read_csv('water_solubility.txt', sep='\t', header=None) # データの読み込み
df.columns = ['SMILES', 'CAS', 'LogS'] #カラム名の変更

#smiles→molファイルの変換
mols_list = []
for smile in df['SMILES']:
    try:
        mols_list.append(Chem.MolFromSmiles(smile))
    except:
        mols_list.append(np.nan)

#molファイルに変換できなかったSMILESの行番号を確認しdfから削除
drop_index = pd.Series(mols_list)[pd.Series(mols_list).isnull()].index.to_list()
df = df.drop(index = drop_index).reset_index(drop = True)

#再度smiles→molファイルの変換
mols_list = []
for smile in df['SMILES']:
    try:
        mols_list.append(Chem.MolFromSmiles(smile))
    except:
        mols_list.append(np.nan)

MOLファイルへの変換まで完了したので、mordred記述子を使って「MOLファイル→ベクトル」へ変換していきます。

「MOLファイル→ベクトル」への変換

mordred記述子を使って「MOLファイル→ベクトル」へ変換します。
下記のコードでは、mols_listの中にあるMOLファイルをベクトルに変換し、Pandas DataFrameの形で出力し、X_trainに出力しています。
X_trainを出力すると「MOLファイル→1826次元のベクトル」に変換されていることが確認できます。

from mordred import Calculator, descriptors
calc = Calculator(descriptors)
X_train = calc.pandas(mols_list)

ベクトルに変換できたので、続いて欠損値処理処理方法について解説します。

mordred記述子で変換後のX_trainの欠損値処理

mordred記述子を使って「MOLファイル→ベクトルに変換」した場合、mordred記述子で変換できなかった箇所（欠損箇所）は、英語（str型）で表記されることに注意してください。
X_trainをいつものようにX_train.isnull().sum()としても、欠損値は0（欠損値はない）と表示されますが、欠損値処理が必要です。
X_trainの英語で記載されている箇所は欠損値です。

欠損値処理のやり方としては、欠損箇所（英語の箇所）を一旦欠損値（np.nan）に変換し、その後欠損値を処理します。まずは欠損箇所（英語の箇所）を一旦欠損値（np.nan)に変換します。
下記のコードでは、X_trainを’str’型に変換し、文字’a’~’z’のいずれかを含むデータを欠損値（np.nan）に変換し、その後float型に変換しています。
こうすることで英語を含むデータは全て欠損値に変換できます。
（この書き方が最善かはわからないですが、一応これで処理できます。）

X_train = X_train.astype('str')
X_train = X_train[~X_train.apply(lambda d: d.str.contains('a|b|c|d|e|f|g|h|i|j|k|l|m|n|o|p|q|r|s|t|u|v|w|x|y|z',na = True))]
X_train = X_train.astype(float)
X_train.isnull().sum()[X_train.isnull().sum() > 0]

欠損箇所（英語の箇所）を欠損値(np.nan)に変換できたので、欠損値が多い列を削除します。
今回は列のデータのうち5%以上欠損値がある列を削除してみます。
列数が1826→1195となりました。
（５％が正解というわけではないので、色々試してください。）

defect_rate = X_train.isnull().sum()[X_train.isnull().sum() >=0] / len(X_train)
columns = defect_rate[defect_rate < 0.05].index.tolist()
X_train = X_train[columns]
X_train

最後に、欠損値をその列の平均値で補完します。
欠損値補完後、欠損値の数が0であることが確認できます。
（欠損値の補完方法は列の平均値以外にも色々考えられるので、遊びで色々試してください。予測精度が変わると思います）

for column in X_train.columns:
    mean = X_train[column].mean()
    X_train[column] = X_train[column].fillna(mean)
    
X_train.isnull().sum()[X_train.isnull().sum() > 0]

X_trainの欠損値処理が完成（mordred記述子を使って化学構造をベクトルに変換）したので、モデル構築していきます。

モデル構築

計算のため、上で算出したX_trainとy_trainを結合してdf_trainとします。

y_train = df['LogS']
df_train = pd.concat([X_train,y_train], axis = 1)

それではデータの前処理が終わったので学習を実施し、予測精度の確認をしていきます！

学習と予測精度の確認

今回もモデルはLightGBMを使用しています。
精度検証はKFoldを使って10分割交差検証を実施しています。結果は、y_test（正解データ）とy_pred(モデルの予測結果）をverificationに格納しています。
verificationを出力すると、y_testとy_predが格納されていることがわかります。

from sklearn.model_selection import KFold
import lightgbm as lgb

verification = pd.DataFrame()
verification['y_test'] = df['LogS']

#KFoldの設定
kf = KFold(n_splits = 10, shuffle = True) 
#交差検証
for train_idx, test_idx in kf.split(df):
    X_train_ = df_train.iloc[train_idx, :-1]
    y_train_ = df_train.iloc[train_idx, -1]
    X_test = df_train.iloc[test_idx, :-1]
    y_test = df_train.iloc[test_idx, -1]
    
    model = lgb.LGBMRegressor()
    model.fit(X_train_, y_train_)
    verification.loc[test_idx, 'y_pred'] = model.predict(X_test)

verification

最後に予測精度結果を下記に示します。
グラフから全てのyの値で高い予測精度であることが確認できます。
またr2_scoreは0.92と非常に高い予測精度である事が確認できます。
ハイパーパラメータを調整することで、もう少し予測精度を上げれそうです。

import matplotlib.pyplot as plt
plt.figure(figsize=(5,5))

plt.scatter(verification['y_test'], verification['y_pred']) # 散布図を描画

plt.show()

from sklearn.metrics import r2_score
r2_score(verification['y_test'], verification['y_pred'])