【ケモインフォマティクス入門：コード付き】化合物の溶解度予測①【データ取得〜記述子作成】
今回はケモインフォマティクス入門編ということで、化合物の溶解度予測を行います。
本当は各ステップで説明したいことがありますが、読みやすさを考えて、新しい記事を作成してリンクを貼る方式で進めようと考えています。
細かい説明を済ませて足場を固めてから今回のような実践的な内容に移っても良いですが、「とりあえずやってみる」というのが初見の分野ではおすすめです。
コード付きで説明していきますので、あまりプログラミングの経験がない方でも試せる内容になっています。
ぜひお試しください！
溶解度予測の大まかな流れ
今回、化合物の溶解度予測をする流れは以下の通りです。
1. データセットの取得
2. RDKitを用いてSMILESを記述子に変換
3. データの前処理
4. Scikit-learnを用いて記述子と溶解度の関係を学習
5. 学習に用いていないデータに対する溶解度の予測精度の評価
本記事では1~2. について説明します。
1. データセットの取得
溶解度予測モデルの学習や予測精度の評価に用いるためのデータを取得します。
実務で用いる場合は、社内に存在するデータや論文・特許などの文献値データをExcelなどで収集することが多いかと思います。
ここではそういったデータに対する包括的な説明はできないので、ネット上で公開されているデータを採用します。
今回は下記の論文に付属するデータセットを用います。
「ESOL: Estimating Aqueous Solubility Directly from Molecular Structure」
J. Chem. Inf. Comput. Sci., 2004, 44, 1000-1005.
Supporinf Information Availableの中に化合物の溶解度データセットがコンマ区切りデータとしてテキストファイルで公開されています。
csvファイルで保存するならば、テキストファイルを開いてNumbersなどに貼り付けてcsv形式で保存すればできますね。
どのようなデータセットなのかを確認してみましょう。
import pandas as pd df = pd.read_csv('solubility_dataset.csv') df = pd.read_csv('solubility_dataset.csv') df.head()
サンプル数は1144個の化合物があり、それぞれに対して以下の情報が与えられていることがわかります。
- Compound ID
- measured log(solubility:mol/L) → 溶解度測定値
- ESOL predicted log(solubility:mol/L) → 論文中での溶解度予測値
- SMILES
2. RDKitを用いてSMILESを記述子に変換
SMILESは文字列の形式なので、データ予測には不向きです。コンピュータの扱いやすい記述子(数字)の形式に変換しましょう。
今回はRDKitという化学分野で機械学習の内容を扱うソフトウェアのようなものを使ってSMILESを変換します。RDKitの処理により記述子を200個生成します。
この部分の処理はネット上で検索してもあまり丁寧に説明している記事が無い印象です。RDKitの処理部分んフォーカスした記事もいずれ作っていこうと思います。
from rdkit import Chem from rdkit.Chem import AllChem, Draw, PandasTools, Descriptors from rdkit.ML.Descriptors import MoleculeDescriptors #RDKitでSMILES読み込んで、sdf作成 IDs = df["Compound ID"] measured_solubilities = df["measured log(solubility:mol/L)"] predicted_solubilities = df["ESOL predicted log(solubility:mol/L)"] SMILES = df["SMILES"] mols = [Chem.MolFromSmiles(SMILES) for SMILES in SMILES] for mol in mols: Chem.AddHs(mol) AllChem.Compute2DCoords(mol) writer = Chem.SDWriter('solubility_dataset.sdf') for (mol,ID, measured_solubility, predicted_solubility ,SMILES) in zip(mols,IDs,measured_solubilities, predicted_solubilities, SMILES): mol.SetProp("ID",str(ID)) mol.SetProp("measured_solubility",str(measured_solubility)) mol.SetProp("predicted_solubility",str(predicted_solubility)) mol.SetProp("SMILES",str(SMILES)) writer.write(mol) writer.close() #RDKitでsdf読み込んで、記述子に変換 mols = Chem.SDMolSupplier('solubility_dataset.sdf') IDs = "Compound ID" measured_solubilities = "measured log(solubility:mol/L)" predicted_solubilities = "ESOL predicted log(solubility:mol/L)" SMILES = "SMILES" mols = [mol for mol in mols if mol is not None] #入力構造に不備があるとNoneになるため descriptor_names = [descriptor_name[0] for descriptor_name in Descriptors._descList] desc_calc = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names) descriptors = pd.DataFrame([desc_calc.CalcDescriptors(mol) for mol in mols]) descriptors.columns = descriptor_names #元のdfにくっつけて、SMILES削除、Nameをindex df = pd.concat([df, descriptors], axis=1) df = df.drop(['SMILES'], axis=1) df = df.set_index("Compound ID") df.to_csv('solubility_dataset_Descriptor.csv') df.head()
上記のようなコードでSMILESからsdfに変換し、sdfから記述子に変換して元のデータにくっつけるような操作ができます。dataframeの列数が202となっており、最初の2つは実測溶解度、計算溶解度なので、SMILESから生成した記述子は200個になることが分かります。
RDKitの使い方に慣れていないとややおまじない的な感じになってしまいますね。この記事に肉付けする形で説明を施していきたいと思います。
最後に説明変数と目的変数を分離して本記事をまとめようと思います。
features = df.iloc[:, 2:] target = df['ESOL predicted log(solubility:mol/L)'] print('shape of features: {}'.format(features.shape)) print('shape of target: {}'.format(target.shape))
上記の操作で説明変数は1144化合物に対して200の特徴量を持ち、目的変数は実測の溶解度とすることができました。
まとめ
本記事では溶解度データセットをダウンロードし、SMILESを200個の説明変数に変換する処理を行いました。
今回は「まずは試してみる」を目標にしたので、各段階で説明しなければよくわからない部分もあったかと思います。
今後の記事でそういった下準備的な知識も補完できるようにしていきたいと思いますので、どうぞよろしくお願いします。
