目的

化学を専門とする機械学習の初学者が書くケモインフォマティクスの記事です。

pythonを用いてサポートベクターマシンによるごく簡単な機械学習を実装します。
有機化学分野における機械学習の概観を掴むことが目標です。

化学における機械学習

近年、他のあらゆる分野と同様に、化学分野においても機械学習の活用が活発に行われています。機械学習は、教師あり学習、教師なし学習、強化学習に分類されます。教師あり学習は画像解析や文字認識、教師なし学習はクラスタリングやマーケティング、強化学習はゲームAIや自動運転などに応用されます。化学分野の実用例は、教師あり学習はデータベースを用いた物性予測、教師なし学習はスペクトル解析、強化学習は実験・合成条件の最適化などがあり、実際に論文として報告されています。

本記事では、化学における機械学習を極々単純化した例として、分子が芳香族であるかどうかを判定するモデルをpythonで実装します。このモデルは教師あり学習であり、サポートベクターマシンというアリゴリズムで学習を行います。この機械学習の流れは以下の通りです。

I. 教師データの収集
II. データの前処理
III. 学習
IV. 予測

機械学習の実装

分子が芳香族であるかどうかを判定するモデルを実装します。
Google Colaboratoryで動作することを確認しています。
今回使用するモジュールを先にimportしておきます。

!pip install rdkit
!pip install mordred
!pip install numpy==1.23.5

from rdkit import Chem
from rdkit.Chem import AllChem, Draw, rdMolDescriptors, MACCSkeys
from rdkit.Chem.Descriptors import MolWt, MolLogP
from rdkit.Avalon.pyAvalonTools import GetAvalonFP
from mordred import Calculator, descriptors
import pandas as pd

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
from sklearn.metrics import accuracy_score, f1_score

from mlxtend.plotting import plot_decision_regions
import numpy as np

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

I. 教師データの収集

学習をするにあたって、学習の元となる教師データを準備します。今回はごく簡単な例なので以下の6つの分子を教師データとします。

これらの分子は芳香族かどうかをラベリングをしておきます。

2-メチルブタン、芳香族でない
プロピレン、芳香族でない
アセトン、芳香族でない
ベンゼン、芳香族である
m-キシレン、芳香族である
安息香酸、芳香族である

最終的に知りたい芳香族かどうかを目的変数と呼びます。これをyとして定義しておきます。

y = np.array([0, 0, 0, 1, 1, 1])  # 目的変数

II. データの前処理

ここで分子構造を機械学習で扱いやすいデータに変換します。手順としては、以下となります。

i. 分子構造をSMILESに書き換え
ii. SMILESから分子記述子を生成

i. 分子構造をSMILESに書き換え

まずSMILESはsimplified molecular input line entry systemの略称で、分子を文字列化したものです。文字列として計算機に入力できるのでケモインフォマティクスの分野でよく用いられています。以下の手順で分子をSMILESに書き換えることができます。

SMILESの書き方

a. 水素原子(H)を省略する。
b. 単結合を省略し、二重結合を” = “、三重結合を” # ”で表す。
c. 分岐は丸括弧で表現する。
d. 環状分子は環を巻く結合をもつ原子の後ろに同一番号をつける。芳香環を構成する元素は小文字にする。

このルールに則って教師データをSMILESに書き換えていきます。

smiles_list = [
    'CC(C)CC',  # 2-メチルブタン
    'CC=C',  # プロピレン
    'CC(=O)C',  # アセトン
    'c1ccccc1',    # ベンゼン
    'Cc1cccc(C)c1',  # m-キシレン
    'c1ccc(c(c1)C(=O)O)',  # 安息香酸
]

ii. SMILESから分子記述子を生成

分子記述子とは、分子をプログラムで扱えるようにした構造や物性などの特徴量のベクトルデータのことです。多種多様な分子記述子が考案されています。例えば、先ほどの教師データをMordred記述子に変換すると以下のようになります。

Mordred記述子は分子構造とそれから導かれる計算値を1613要素の配列にしたものです。今回はごく簡単な例として、以下の分子記述子をXとして用いてみます。

X = [分子量、水素の数]

Mordred記述子に両者の値が含まれていますのでそれを利用します。このような、学習に用いる変数を説明変数と呼びます。

# 分子オブジェクトを生成
molecules = [Chem.MolFromSmiles(smile) for smile in smiles_list]

# Mordred記述子
calc = Calculator(descriptors, ignore_3D=True)
mordred_descriptors = calc.pandas(molecules)

# Mordred記述子をデータフレームに変換
mordred_df = mordred_descriptors
mordred_df.index = smiles_list

X = mordred_df[['MW', 'nH']].values  # 説明変数

III. 学習

学習アルゴリズムにsupport vector machine (SVM)を用いて学習を行っていきます。SVMは二つのクラスのデータ群を分割するような境界を決定するアルゴリズムです。わかりやすいようにそれを図示します。

先にSVMの結果を図示する関数を定義します。

# 決定領域表示関数
def plot_decision_regions(X, y, classifier, test_idx=None, resolution=0.02):

    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    x1_min, x1_max = X[:, 0].min() - 15, X[:, 0].max() + 15
    x2_min, x2_max = X[:, 1].min() - 3, X[:, 1].max() + 3
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.3, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # 全てのサンプルを表示
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0],
                    y=X[y == cl, 1],
                    alpha=0.8,
                    c=colors[idx],
                    marker=markers[idx],
                    label=cl,
                    edgecolor='black')

学習は次のように実行します。

# 目的変数の配列を2次元に変換
y = y.reshape(len(y), 1)

# 目的変数をstr型→int型に変換
label_names = list(dict.fromkeys(y[:, 0]))
label_dict = dict(zip(label_names, range(len(label_names))))
y_int=np.vectorize(lambda x: label_dict[x])(y)

# 学習
clf = Pipeline([
    ("scaler", StandardScaler()),
    ("liner_svc", LinearSVC(C=1, loss="hinge")),
])
clf.fit(X, y_int)

# mlxtendで決定境界可視化
plot_decision_regions(X, y_int[:, 0], classifier=clf, test_idx=None, resolution=0.02)

plt.xlabel('Molecular Weight')
plt.ylabel('Number of H')

学習結果を図示したものが以下になります。

図中の番号が先に示した分子の番号を示しています。今回は分子記述子をXが単純なので確認してみますと、1番の2-メチルブタンは分子量が72 g/mol、水素の数が12個なので図の上の方に来ていて、4番のベンゼンは分子量が78 g/mol、水素の数が6個なので下の方に来ています。良さそうです。

説明変数が2種類の時のSVMは二次元平面に領域を分割する直線を引くことと対応します。この直線$g(x,y)$がモデルであり、未知のデータに対しては$g(x,y)$との大小比較により目的変数の予測を行います。図中には、説明変数が0の芳香族でないデータが赤四角、説明変数が1の芳香族であるデータが青バツで示されています。今回、説明変数がたまたま良かったので綺麗にモデルを導出することができました。

実際のSVMによる機械学習では説明変数のn次元の空間を分割する境界がモデルとなります。

IV. 予測

最後に予測を行います。本当は芳香族かどうかが未知の分子が良いですが、今回はごく簡単な例なので、芳香族ということがわかっているナフタレンについて予測しましょう。予測結果として1が出て来れば正解です。

では、ナフタレンについて、教師データと同様のデータの前処理を行います。

smiles_list2 = [
    'c1ccc2c(c1)cccc2',  #ナフタレン
]

# 分子オブジェクトを生成
molecules2 = [Chem.MolFromSmiles(smile) for smile in smiles_list2]

# Mordred記述子
calc2 = Calculator(descriptors, ignore_3D=True)
mordred_descriptors2 = calc2.pandas(molecules2)

# Mordred記述子をデータフレームに変換
mordred_df2 = mordred_descriptors2
mordred_df2.index = smiles_list2

X2 = mordred_df2[['MW', 'nH']].values  # 説明変数

予測は次のように行います。

y_pred = clf.predict(X2)

結果は次のようになります

array([1])

今回のモデルはナフタレンに対しては正解を出しました。
実際にプロットしてみるとモデルの直線の下にくるので、うまく予測ができたことがわかります。

まとめ

pythonを用いてSVMによるごく簡単な有機化学に関する機械学習を実装しました。

今回は、結果的にナフタレンが芳香族であることをうまく予測することができましたが、これは今回用意した説明変数Xに対して教師データが適切だったからです。例えば、不飽和度の大きな芳香族でない炭化水素が教師データに含まれていたら、芳香族でないのにも関わらず学習結果の図の下の方に来るので、説明変数Xのみではうまく学習できていなかったかもしれません。

実際の学習では教師データの数を増やし、説明変数の次元を増やし、種々のアルゴリズムを活用して精度の高い機械学習を行います。今回の学習は実用レベルにはないですが、概観を掴むにはわかりやすいかと思います。