化学情報学にRDKitとScikit-learnで入門する


このページを見たあなたは、もう化学情報学（ケモインフォマティクス）という世界に足を踏み入れている。まずはこの世界の入り口を概観してみよう。

# 題材として、水溶解度を取り扱ってみます

Water solublity (logS) database http://modem.ucsd.edu/adme/databases/databases_logS.htm にある「database.dat (the SMILES file includes 1290 molecules)」 を使ってみましょう。

### まず、データファイルの URL をチェックします。

http://modem.ucsd.edu/adme/data/databases/logS/data_set.dat であることが分かります。今回の目的は、このファイルの左端の列にある文字列（SMILES）から、水溶解度（logS）を予測する機械学習モデルを作成することとします。

### データファイルのダウンロード

手でクリクリしてダウンロードしてもいいですが、せっかくなので urllib.request を使ってみましょう。


```python
import urllib.request # URLで指定したファイルをダウンロードするライブラリ

url = 'http://modem.ucsd.edu/adme/data/databases/logS/data_set.dat'
urllib.request.urlretrieve(url, 'water_solubility.txt') 
```




    ('water_solubility.txt', <http.client.HTTPMessage at 0x109eb8438>)



ファイルを water_solubility.txt という名前で保存しました。

### ダウンロードしたデータファイルを pandas で読み込む

読み込み方の詳細は <a href="https://qiita.com/maskot1977/items/1f6c51965b1bfc341e4b">タブ区切りデータ、コンマ区切りデータ等の読み込み</a> をご参照ください。


```python
import pandas as pd

df1 = pd.read_csv('water_solubility.txt', sep='\t', header=None) # データの読み込み
df1 # 読み込んだデータの確認
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>O=C(C)N</td>
      <td>60-35-5</td>
      <td>1.58</td>
    </tr>
    <tr>
      <th>1</th>
      <td>NNC</td>
      <td>60-34-4</td>
      <td>1.34</td>
    </tr>
    <tr>
      <th>2</th>
      <td>O=C(C)O</td>
      <td>64-19-7</td>
      <td>1.22</td>
    </tr>
    <tr>
      <th>3</th>
      <td>N1CCCC1</td>
      <td>123-75-1</td>
      <td>1.15</td>
    </tr>
    <tr>
      <th>4</th>
      <td>O=C(N)NO</td>
      <td>127-07-1</td>
      <td>1.12</td>
    </tr>
    <tr>
      <th>5</th>
      <td>OCC(O)CO</td>
      <td>56-81-5</td>
      <td>1.12</td>
    </tr>
    <tr>
      <th>6</th>
      <td>O=C(C)N(C)C</td>
      <td>127-19-5</td>
      <td>1.11</td>
    </tr>
    <tr>
      <th>7</th>
      <td>[nH0]1cccc[nH0]1</td>
      <td>289-80-5</td>
      <td>1.10</td>
    </tr>
    <tr>
      <th>8</th>
      <td>[nH0]1ccc[nH0]c1</td>
      <td>289-95-2</td>
      <td>1.10</td>
    </tr>
    <tr>
      <th>9</th>
      <td>OCC(C(C(C(CO)O)O)O)O</td>
      <td>50-70-4</td>
      <td>1.09</td>
    </tr>
    <tr>
      <th>10</th>
      <td>O=C(CC(C)N)O</td>
      <td>541-48-0</td>
      <td>1.08</td>
    </tr>
    <tr>
      <th>11</th>
      <td>N1CCNCC1</td>
      <td>110-85-0</td>
      <td>1.07</td>
    </tr>
    <tr>
      <th>12</th>
      <td>N1CCCCC1</td>
      <td>110-89-4</td>
      <td>1.07</td>
    </tr>
    <tr>
      <th>13</th>
      <td>c1c(O)cc[nH0]c1</td>
      <td>626-64-2</td>
      <td>1.02</td>
    </tr>
    <tr>
      <th>14</th>
      <td>c1c(O)[nH0]ccc1</td>
      <td>72762-00-6</td>
      <td>1.02</td>
    </tr>
    <tr>
      <th>15</th>
      <td>O=C(O)CCCC(O)=O</td>
      <td>110-94-1</td>
      <td>1.00</td>
    </tr>
    <tr>
      <th>16</th>
      <td>C1N(C)CCOC1</td>
      <td>109-02-4</td>
      <td>1.00</td>
    </tr>
    <tr>
      <th>17</th>
      <td>O=C(N)OC</td>
      <td>598-55-0</td>
      <td>0.97</td>
    </tr>
    <tr>
      <th>18</th>
      <td>CN1C(=O)C=CC=C1</td>
      <td>694-85-9</td>
      <td>0.96</td>
    </tr>
    <tr>
      <th>19</th>
      <td>O=C(O)CCl</td>
      <td>79-11-8</td>
      <td>0.93</td>
    </tr>
    <tr>
      <th>20</th>
      <td>O=C(N)OCC</td>
      <td>51-79-6</td>
      <td>0.85</td>
    </tr>
    <tr>
      <th>21</th>
      <td>COC(=O)C1=CCCN(C)C1</td>
      <td>63-75-2</td>
      <td>0.81</td>
    </tr>
    <tr>
      <th>22</th>
      <td>c1c(O)cc(cc1)O</td>
      <td>108-46-3</td>
      <td>0.81</td>
    </tr>
    <tr>
      <th>23</th>
      <td>c1cc[nH0]cc1</td>
      <td>110-86-1</td>
      <td>0.76</td>
    </tr>
    <tr>
      <th>24</th>
      <td>N1CCNCC1C</td>
      <td>109-07-9</td>
      <td>0.74</td>
    </tr>
    <tr>
      <th>25</th>
      <td>OCC1(C(C(C(CO1)O)O)O)O</td>
      <td>@57-48-7</td>
      <td>0.64</td>
    </tr>
    <tr>
      <th>26</th>
      <td>c1c(O)c(O)ccc1</td>
      <td>120-80-9</td>
      <td>0.62</td>
    </tr>
    <tr>
      <th>27</th>
      <td>O=C(N)c1ccc[nH0]c1</td>
      <td>98-92-0</td>
      <td>0.61</td>
    </tr>
    <tr>
      <th>28</th>
      <td>C=1C(=O)NC=NC1</td>
      <td>51953-17-4</td>
      <td>0.59</td>
    </tr>
    <tr>
      <th>29</th>
      <td>O=CCC</td>
      <td>123-38-6</td>
      <td>0.58</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>1260</th>
      <td>c1c(Cl)ccc(c1)c1ccc(cc1Cl)Cl</td>
      <td>7012-37-5</td>
      <td>-6.21</td>
    </tr>
    <tr>
      <th>1261</th>
      <td>ClC12C3C(C4C(C3Cl)O4)C(C(=C1Cl)Cl)(C2(Cl)Cl)Cl</td>
      <td>1024-57-3</td>
      <td>-6.29</td>
    </tr>
    <tr>
      <th>1262</th>
      <td>O=C(C1C(C=C(Cl)Cl)C1(C)C)OCc1cccc(c1)Oc1ccccc1</td>
      <td>52645-53-1</td>
      <td>-6.29</td>
    </tr>
    <tr>
      <th>1263</th>
      <td>Clc1c(Cl)ccc(c1c1ccccc1)Cl</td>
      <td>55702-45-9</td>
      <td>-6.29</td>
    </tr>
    <tr>
      <th>1264</th>
      <td>ClC1(C2(C3C=CC(C3C1(C(=C2Cl)Cl)Cl)Cl)Cl)Cl</td>
      <td>76-44-8</td>
      <td>-6.32</td>
    </tr>
    <tr>
      <th>1265</th>
      <td>Clc1c(c2c(Cl)cc(cc2)Cl)ccc(c1)Cl</td>
      <td>2437-79-8</td>
      <td>-6.51</td>
    </tr>
    <tr>
      <th>1266</th>
      <td>Clc1ccccc1C(c1ccc(cc1)Cl)C(Cl)(Cl)Cl</td>
      <td>789-02-6</td>
      <td>-6.62</td>
    </tr>
    <tr>
      <th>1267</th>
      <td>Cc1cccc(c1O[P+]([O-])(Oc1ccc(cc1C)O)Oc1cc(ccc1...</td>
      <td>1330-78-5</td>
      <td>-6.70</td>
    </tr>
    <tr>
      <th>1268</th>
      <td>O=C(CCCCCCCCCCCCCCC)O</td>
      <td>57-10-3</td>
      <td>-6.81</td>
    </tr>
    <tr>
      <th>1269</th>
      <td>c1c(C(C(Cl)(Cl)Cl)c2ccc(cc2)OC)ccc(c1)OC</td>
      <td>72-43-5</td>
      <td>-6.89</td>
    </tr>
    <tr>
      <th>1270</th>
      <td>O=C(c1ccccc1C(=O)OCC(CC)CCCC)OCC(CC)CCCC</td>
      <td>117-81-7</td>
      <td>-6.96</td>
    </tr>
    <tr>
      <th>1271</th>
      <td>Clc1c(c2c(Cl)ccc(c2)Cl)cc(cc1)Cl</td>
      <td>35693-99-3</td>
      <td>-7.00</td>
    </tr>
    <tr>
      <th>1272</th>
      <td>c1c(c2c(Cl)c(c(cc2)Cl)Cl)c(Cl)c(cc1)Cl</td>
      <td>52663-62-4</td>
      <td>-7.05</td>
    </tr>
    <tr>
      <th>1273</th>
      <td>c1c(c2c(Cl)c(Cl)c(c(c2)Cl)Cl)cccc1</td>
      <td>33284-53-6</td>
      <td>-7.16</td>
    </tr>
    <tr>
      <th>1274</th>
      <td>OCCCCCCCCCCCCCCCC</td>
      <td>36653-82-4</td>
      <td>-7.26</td>
    </tr>
    <tr>
      <th>1275</th>
      <td>ClC1=C(C2(C3C4CC(C=C4)C3C1(C2(Cl)Cl)Cl)Cl)Cl</td>
      <td>309-00-2</td>
      <td>-7.33</td>
    </tr>
    <tr>
      <th>1276</th>
      <td>c1c2c3c4c(cccc4c4c2cccc4)ccc3cc1</td>
      <td>192-97-2</td>
      <td>-7.80</td>
    </tr>
    <tr>
      <th>1277</th>
      <td>c1c(Cl)ccc(c1c1ccc(c(Cl)c1Cl)Cl)Cl</td>
      <td>38380-02-8</td>
      <td>-7.91</td>
    </tr>
    <tr>
      <th>1278</th>
      <td>c1c2Cc3c(c2ccc1)cc1c(cccc1)c3</td>
      <td>243-17-4</td>
      <td>-8.04</td>
    </tr>
    <tr>
      <th>1279</th>
      <td>c1c2c(ccc3c4c(ccc23)cccc4)ccc1</td>
      <td>218-01-9</td>
      <td>-8.06</td>
    </tr>
    <tr>
      <th>1280</th>
      <td>Clc1c(c2c(Cl)c(Cl)c(cc2)Cl)c(Cl)c(c(c1)Cl)Cl</td>
      <td>52663-71-5</td>
      <td>-8.30</td>
    </tr>
    <tr>
      <th>1281</th>
      <td>c1c(Cl)c(Cl)cc(c1c1ccc(c(Cl)c1Cl)Cl)Cl</td>
      <td>35065-28-2</td>
      <td>-8.32</td>
    </tr>
    <tr>
      <th>1282</th>
      <td>c1c(Cl)c(Cl)c(c(c1c1cccc(Cl)c1Cl)Cl)Cl</td>
      <td>55215-18-4</td>
      <td>-8.42</td>
    </tr>
    <tr>
      <th>1283</th>
      <td>c1c2c3c(cccc3c3c2cc2c(c3)cccc2)cc1</td>
      <td>207-08-9</td>
      <td>-8.49</td>
    </tr>
    <tr>
      <th>1284</th>
      <td>c1c2cc3c(cc4c(c3)cccc4)cc2ccc1</td>
      <td>92-24-0</td>
      <td>-8.60</td>
    </tr>
    <tr>
      <th>1285</th>
      <td>Clc1c(c2c(ccc(Cl)c2Cl)Cl)c(ccc1Cl)Cl</td>
      <td>38411-22-2</td>
      <td>-8.65</td>
    </tr>
    <tr>
      <th>1286</th>
      <td>c1c2c(ccc3c2cccc3)cc2c1ccc1ccccc12</td>
      <td>53-70-3</td>
      <td>-8.66</td>
    </tr>
    <tr>
      <th>1287</th>
      <td>c1c2c3cccc4ccc5ccc6c(c2c(cc6)cc1)c5c43</td>
      <td>191-24-2</td>
      <td>-9.03</td>
    </tr>
    <tr>
      <th>1288</th>
      <td>Clc1c(Cl)c(c2cc(c(c(Cl)c2Cl)Cl)Cl)c(c(c1Cl)Cl)Cl</td>
      <td>40186-72-9</td>
      <td>-10.26</td>
    </tr>
    <tr>
      <th>1289</th>
      <td>Clc1c(c2c(Cl)c(Cl)c(c(c2Cl)Cl)Cl)c(Cl)c(cc1Cl)Cl</td>
      <td>52663-77-1</td>
      <td>-10.41</td>
    </tr>
  </tbody>
</table>
<p>1290 rows × 3 columns</p>
</div>



このようにして、1290行3列の行列が得られました。今回の目的は、この「0列目」の文字列（SMILES）から「2列目」の数字（logS）を予測する問題を解く、ということになります。

# RDKit を使ってみる

まずは化合物の絵を描いてみたい。


```python
from rdkit import Chem
```


```python
mol = Chem.MolFromSmiles("O=C(N)c1ccc[nH0]c1")
```


```python
from rdkit.Chem import Draw
Draw.MolToImage(mol)
```


![output_7_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/56b11cfe-1749-e2bb-5ca0-5890f5abc0df.png)



化合物の絵が描けることを確認しました。今度は、今回ダウンロードした水溶解度のデータに含まれている化合物の絵を描きましょう。ただし、全部描くのはしんどいので、先頭いくつかだけ。


```python
mols = []
for smile in df1[0]:
    mols.append(Chem.MolFromSmiles(smile))
```


```python
# 上と同じ意味。「リスト内包表記」という書き方です。
mols = [Chem.MolFromSmiles(smile) for smile in df1[0]]
```


```python
Draw.MolsToGridImage(mols[100:124], molsPerRow=5, subImgSize=(200,200))
```


![output_11_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/3e1df1de-3d26-ea3c-7dfe-a06b5079af89.png)



## 化合物フィンガープリント

化合物を数字で表現する方法のひとつとして、フィンガープリントと呼ばれるベクトルがあります。今回はその一種であるモーガンフィンガープリントを使ってみましょう。


```python
from rdkit.Chem import AllChem
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048)
```


```python
fingerprint = []
for s in fp.ToBitString():
    fingerprint.append(int(s))
```


```python
# 上と同じ意味
fingerprint = [int(s) for s in fp.ToBitString()]
```

水溶解度のデータに含まれる化合物をフィンガープリントにします。ただし、エラーが出て計算できない化合物がある場合は除外します。


```python
import numpy as np
fingerprints = []
safe = []
for mol_idx, mol in enumerate(mols):
    try:
        fingerprint = [x for x in AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048)]
        fingerprints.append(fingerprint)
        safe.append(mol_idx)
    except:
        print("Error", mol_idx)
        continue
fingerprints = np.array(fingerprints)
```

    Error 175



```python
print(fingerprints.shape)
pd.DataFrame(fingerprints).head()
```

    (1289, 2048)





<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
      <th>...</th>
      <th>2038</th>
      <th>2039</th>
      <th>2040</th>
      <th>2041</th>
      <th>2042</th>
      <th>2043</th>
      <th>2044</th>
      <th>2045</th>
      <th>2046</th>
      <th>2047</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>...</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>...</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>...</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>...</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>...</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
      <td>0</td>
    </tr>
  </tbody>
</table>
<p>5 rows × 2048 columns</p>
</div>



# フィンガープリントから水溶解度を回帰する

回帰する方法はいくつかありますが、ランダムフォレスト というものを使ってみます。


```python
X = fingerprints
Y = df1.iloc[safe, 2]
```


```python
from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor()
rf.fit(X,Y)

from sklearn.metrics import r2_score
print("R2=", r2_score(Y, rf.predict(X)))

import numpy as np
from sklearn.metrics import mean_squared_error
print("RMSE=", np.sqrt(mean_squared_error(Y, rf.predict(X))))

from sklearn.metrics import mean_absolute_error
print("MAE=", mean_absolute_error(Y, rf.predict(X)))

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(Y, rf.predict(X), alpha=0.2, c="blue")
plt.plot([Y.min(), Y.max()], [Y.min(), Y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    /Users/kot/miniconda3/envs/py3new/lib/python3.6/site-packages/sklearn/ensemble/forest.py:248: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.
      "10 in version 0.20 to 100 in 0.22.", FutureWarning)


    R2= 0.9381194414671183
    RMSE= 0.5054545837347694
    MAE= 0.34466386277058253

![output_21_2.png](https://qiita-image-store.s3.amazonaws.com/0/211162/86968b98-3ce2-dcaf-7372-af45437b35c6.png)


うん、あまり混みいったことしてない割には、いい感じ。では、サポートベクトル回帰ってやつも使ってみましょうか。


```python
from sklearn.svm import SVR
svr = SVR()
svr.fit(X,Y)

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(Y, svr.predict(X), alpha=0.2, c="blue")
plt.plot([Y.min(), Y.max()], [Y.min(), Y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    /Users/kot/miniconda3/envs/py3new/lib/python3.6/site-packages/sklearn/svm/base.py:196: FutureWarning: The default value of gamma will change from 'auto' to 'scale' in version 0.22 to account better for unscaled features. Set gamma explicitly to 'auto' or 'scale' to avoid this warning.
      "avoid this warning.", FutureWarning)



![output_23_1.png](https://qiita-image-store.s3.amazonaws.com/0/211162/4325eddf-8135-04b3-69f5-1ff6be026217.png)


ありゃ。だめだこりゃ。サポートベクトルマシンは、パラメーター調整をグリグリ頑張って初めて真価を発揮する手法なので、デフォルトパラメータを使うと良くないことが多いのです。

# さあ、次にやることは？

* 化合物の表現方法として、フィンガープリント以外にもさまざまな「化学構造記述子」があります。フィンガープリントにも、いろんなフィンガープリントがあります。

* 機械学習手法として、ランダムフォレスト やサポートベクトルマシン以外にも色々な手法があります。

* パラメータ調節方法として、グリッドサーチやベイズ最適化などの方法があります。

* 機械学習で注意すべき点として「過学習」という現象があり、それを防ぐための方法論として「交差検定」などがあります。

まだまだ学ぶべきことは多いけれども、ここまで読んだあなたは今、こう言える。

## ケモインフォの「世界」に入門したッ！

ゴゴゴゴゴゴ...