ケモインフォマティクスで学ぶscikit-learn

# はじめに

[ケモインフォマティクスで学ぶMatplotlib](https://qiita.com/yukiya285/items/96be7a07122fa63d2713)に引き続き、リピドミクス（脂質の網羅解析）を題材として、Pythonの代表的なライブラリの一つである「Matplotlib」について解説していきます。
ケモインフォマティクスの実践例を中心に説明していきますので、基本を確認したいという人は以下の記事を読んでからこの記事を読んでみてください。

[製薬企業研究者がscikit-learnについてまとめてみた](https://qiita.com/yukiya285/items/75dd6fc539b8c1231404)


# データセットの準備

scikit-learnは、機械学習用のライブラリです。

ここでは、化合物の物性値から、液体クロマトグラフィー（Liquid Chromatography: LC）における保持時間（Retention Time: RT）を、部分最小二乗（partial least squares: PLS）回帰を用いて予測することを考えてみます。

まず、機械学習に用いるデータセットを作成します。

```python
import pandas as pd


params_fatty_acids = ['Heavy atoms', 'Rotatable Bonds', 'van der Waals Molecular Volume', 'logP', 'Molar Refractivity']

lauric = [14, 10, 231.10, 3.99, 59.48]
myristic = [16, 12, 265.70, 4.77, 68.71]
palmitic = [18, 14, 300.30, 5.55, 77.95]
palmitoleic = [18, 13, 297.66, 5.33, 77.85]
stearic = [20, 16, 334.90, 6.33, 87.18]
oleic = [20, 15, 332.26, 6.11, 87.09]
linoleic = [20, 14, 329.62, 5.88, 86.99]
linolenic = [20, 13, 326.98, 5.66, 86.90]
stearidonic = [20, 12, 324.34, 5.44, 86.81]
arachidic = [22, 18, 369.50, 7.11, 96.42]
bishomo_gamma_linolenic = [22, 15, 361.58, 6.44, 96.13]
arachidonic = [22, 14, 358.94, 6.22, 96.04]
eicosapentaenoic = [22, 13, 356.30, 5.99, 95.95]
behenic = [24, 20, 404.10, 7.89, 105.65]
adrenic = [24, 16, 393.54, 7.00, 105.27]
docosapentaenoic = [24, 15, 390.90, 6.77, 105.18]
docosahexaenoic = [24, 14, 388.26, 6.55, 105.09]

df_fatty_acids = pd.DataFrame([lauric, myristic, palmitic, palmitoleic, stearic, oleic, linoleic, linolenic, stearidonic, arachidic, bishomo_gamma_linolenic, arachidonic, eicosapentaenoic, behenic, adrenic, docosapentaenoic, docosahexaenoic], columns=params_fatty_acids)
df_fatty_acids['Experimental Retention Time (min)'] = [4.53, 7.52, 11.02, 10.59, 14.45, 11.86, 9.76, 8.31, 6.71, 17.52, 11.20, 9.96, 8.27, 20.40, 12.75, 11.52, 9.84]

print(df_fatty_acids)
```

ここでは、説明変数として用いる物性パラメーター名のリストを`params_fatty_acids`としています。
それぞれの物性値は[LIPID MAPS](https://www.lipidmaps.org/)のデータベースに格納されている情報を参照しています。
また、RTは、理化学研究所の[PRIMe](http://prime.psc.riken.jp/)というウェブサイトで公開されている[逆相LCにおけるRTデータ](http://prime.psc.riken.jp/Metabolomics_Software/MrmDatabase/Detail%20of%20LCQqQMS%20method%20(ODS-lipids).xlsx)を参照しています。
また、実際は`pandas.read_csv`などでCSVファイルなどを読み込むことが多いと思います。
さらに、欠損値補完などのデータ前処理が必要になることが多いです。


# モデルの構築

次に、予測モデルの構築と、モデルを用いた予測値の算出を行います。

```python
from sklearn.cross_decomposition import PLSRegression


X = df_fatty_acids[params_fatty_acids] # 説明変数
y = df_fatty_acids['Experimental Retention Time (min)'] # 目的変数

pls_rt = PLSRegression()
pls_rt.fit(X, y) # PLS予測モデルを構築

y_pred = pls_rt.predict(X) # 予測値を算出

df_fatty_acids['Predicted Retention Time (min)'] = y_pred
df_fatty_acids['Diff (min)'] = df_fatty_acids['Predicted Retention Time (min)'] - df_fatty_acids['Experimental Retention Time (min)']
df_fatty_acids['Accuracy (%)'] = (df_fatty_acids['Diff (min)'] / df_fatty_acids['Experimental Retention Time (min)']) * 100

print(df_fatty_acids)
```

実測値と予測値の関係を図示してみると、以下のようになります。

```python
%matplotlib inline
import matplotlib.pyplot as plt


plt.scatter(y, y_pred)
plt.xlabel('Experimental Retention Time (min)')
plt.ylabel('Predicted Retention Time (min)')

plt.savefig('rts_fatty_acids.png')
plt.show()
```

![rts_fatty_acids.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/517951/e74072d7-4215-f56a-e65a-9a0317a74992.png)

今回のデータでは実測値と予測値がよく合っていそうですね。
構築したモデルの当てはまり度合いは`r2_score`で確認することができます。

```python
from sklearn.metrics import r2_score


print(r2_score(y, y_pred))
```

`r2_score`は、0から1の間の値をとり、1に近いほど実測値と予測値がよく合っていることを示します。
今回のデータでは、`r2_score`が0.98を越える値となっており、かなり当てはまりの良いモデルということになります。

今回は、RTを予測するのに、5種類の物性パラメータを用いましたが、このうちどれが予測に大きく寄与しているのかを見てみます。

```python
print(pls_rt.coef_)
```

この結果から、今回のデータでは、Rotable Bondsに対する係数（の絶対値）が3.44と最も大きくなっており、この物性値がRTの予測に強く寄与していることが分かります。

# モデルを用いた予測

以上で予測モデル構築に用いたデータの予測精度について議論してきましたが、最後に、モデル構築に用いなかったデータがどれほどの精度で予測できるのか見てみます。

```python
lignoceric = [26, 22, 438.70, 8.67, 114.88]
x_lignoceric = pd.DataFrame([lignoceric], columns=params_fatty_acids)
y_pred_lignoceric = pls_rt.predict(x_lignoceric)

y_exp_lignoceric = 22.31 # 実測値

print(y_exp_lignoceric)
print(y_pred_lignoceric)
```

ここでは、lignoceric acid（FA 24:0）のRTを予測してみました。
予測値と実測値の乖離が1.2分程度となっています。
この差を大きいとみるか小さいとみるかは色々な見方があると思いますが、個人的にはやや予測精度が低いと思います。
理由としては、lignoceric acidが、モデル構築に用いたデータセットに含まれるどの脂肪酸分子種よりも疎水性が高い分子種であり、モデル構築の段階でlignoceric acidの物性に近いデータのfittingが行われていなかったことも関係していると考えられます。

以上のような手順でPLS回帰を行うことができます。
ここでは触れませんでしたが、PLS回帰を行うにあたり、潜在変数の数`n_components`も重要になってきます。
今回はデフォルト値の`2`を用いましたが、ここを変えることで予測精度も少しずつ変わってきます。
また別の機会に解説したいと思います。


# まとめ

ここでは、scikit-learnについて、ケモインフォマティクスで使える実践的な知識を中心に解説しました。
もう一度要点をおさらいしておきましょう。

- scikit-learnを用いることで、機械学習を手軽に行うことができる。
- 機械学習はデータの前処理、モデル構築、予測という流れで行う。
- 構築されたモデルの`r2_score`や、予測値と実測値の差などを見ながら進めていく。


# 参考資料・リンク

[プログラミング言語Pythonとは？AIや機械学習に使える？](https://pharma-engineer.com/basic-programming/python-summary/)