はじめに
前回まで,Kaggle の入門として有名な Titanic のデータセットを用いて機械学習で生存予測に取り組んでいました.
①EDAとロジスティック回帰
②特徴量エンジニアリングによるモデル改善
③tree系モデルの検討
次に,自分の興味があるケモインフォマティクス的な要素についても学びたいと思ったので,化合物の水溶解度について機械学習で予測してみようと思います.
目的変数となるのは水への溶解度 (logS = log₁₀S) です(溶解度 S:mol L⁻¹).
logS = 0 なら 1 mol/L(よく溶ける),logS = −5 なら 0.00001 mol/L(ほとんど溶けない)
という意味です.
今回も AI に学習コースを作ってもらいました.
元となる論文はこちらです.
J. S. Delaney, "ESOL: Estimating Aqueous Solubility Directly from Molecular Structure",
J. Chem. Inf. Comput. Sci., 2004, 44, 1000–1005.
https://pubs.acs.org/doi/10.1021/ci034243x
この論文では,2874 個の化合物の実測溶解度データから,
分子構造の記述子(LogP, 分子量, 芳香族原子の割合, 回転可能結合数)を用いた
線形回帰モデルで溶解度を予測しています.
今回はこのデータセットを使い,段階的にモデルを改善していく流れで学びます.
データセットは DeepChem の GitHub リポジトリから取得します.
https://github.com/deepchem/deepchem/blob/master/datasets/delaney-processed.csv
ライブラリをインポート,データの取得
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from rdkit import Chem
from rdkit.Chem import Draw, Descriptors
# Delaney(ESOL)データセットを取得
url = "https://raw.githubusercontent.com/deepchem/deepchem/master/datasets/delaney-processed.csv"
df = pd.read_csv(url)
データの確認
まず,データの形と列名を調べます.
print("データの形", df.shape)
print("列名", df.columns.tolist())

1128化合物が入っており,各列の意味はざっくり以下の通りです.
- Compound ID: 化合物名
- ESOL predicted...: Delaney の論文の式による予測値(使わない)
- Minimum Degree: 最小結合次数
- Molecular Weight: 分子量
- Number of H-Bond Donors: 水素結合ドナー数
- Number of Rings: 環の数
- Number of Rotatable Bonds: 回転可能結合数
- Polar Surface Area: 極性表面積
- measured log solubility in mols per litre: 実測 logS ← これが予測対象(目的変数)
- smiles: SMILES 文字列
次に,目的変数としてlogS(measured log solubility in mols per litre)を設定し,ヒストグラムを作成します.
target_col = "measured log solubility in mols per litre"
plt.figure(figsize=(8, 5))
plt.hist(df[target_col], bins=40, edgecolor="black")
plt.xlabel("measured logS (mol/L)")
plt.title("Distribution of measured log solubility (n=1128)")
plt.show()

-2あたりに山があり,左(低溶解度)に裾がやや伸びた分布をしていますが,概ね正規分布に近いです.
最初の6分子については構造の描画もしてみます.
smiles列の最初の6つをSMILESからMolオブジェクトへ変換してリストへ入れ,
同様にCompoundIDも取得し,それらの構造を表示させます.
mols = [Chem.MolFromSmiles(smi) for smi in df["smiles"].head(6)]
names = df["Compound ID"].head(6).tolist()
img = Draw.MolsToGridImage(mols, molsPerRow=3, subImgSize=(250, 200), legends=names)
img

天然物であるAmigdalinや農薬のFenfuram
電子材料として使われるPiceneやbenzothiazoleなど,さまざまな分子が存在します.
次に,各記述子とlogSの散布図を作成します.
(Titanicなど分類タスクでは目的変数が0,1だけなので見づらく,使いませんでした)
# 既存の記述子列
desc_cols = ["Minimum Degree", "Molecular Weight", "Number of H-Bond Donors",
"Number of Rings", "Number of Rotatable Bonds", "Polar Surface Area"]
# 各記述子と logS の散布図
fig, axes = plt.subplots(2, 3, figsize=(15, 9))
for ax, col in zip(axes.ravel(), desc_cols): #.ravel()はsubplotsで作った2×3の二次元配列を一次元にする
ax.scatter(df[col], df[target_col], alpha=0.3, s=15)
ax.set_xlabel(col)
ax.set_ylabel("logS")
ax.set_title(f"logS vs {col}")
plt.tight_layout()
plt.show()
加えてヒートマップを作成し,各記述子とlogSの相関と,記述子同士の相関を確認します.
# 既存の記述子 + logS でヒートマップ
all_cols = desc_cols + [target_col]
plt.figure(figsize=(10, 8))
sns.heatmap(df[all_cols].corr(), annot=True, cmap="coolwarm", center=0, fmt=".2f")
plt.title("Correlation between descriptors and logS")
plt.tight_layout()
plt.show()
- 散布図からは,Molecular Weight と Number of Rings が右肩下がりの傾向を示しており,logS との相関が強そうに見えます.
- この二つはヒートマップでもlogSとの相関が-0.64, -0.51と大きな数値なので,中心となる特徴量になりそうです.分子量が大きいほど溶けにくい,環が多いほど溶けにくい,と理に適った説明もできます.
- ヒートマップの記述子同士の相関を見ると,H-Bond Donors と Polar Surface Area の相関が 0.76,Molecular Weight と Number of Rings も 0.65.H-Bond donorになると,極性表面積が増える.大きな分子ほど環が多い,という化学的にも納得できる相関です.
モデル作成・評価
まずは,元論文と同じようにこの6記述子だけでモデルを作成します.線形回帰(Linear regression)を使います.
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import r2_score, mean_squared_error
# 既存の6記述子で予測
desc_cols = ["Minimum Degree", "Molecular Weight", "Number of H-Bond Donors","Number of Rings", "Number of Rotatable Bonds", "Polar Surface Area"]
# X,yの作成
X = df[desc_cols]
y = df[target_col]
# 訓練/テスト分割
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
# パイプライン(標準化 + 線形回帰)
pipe = Pipeline([
("scaler", StandardScaler()),
("model", LinearRegression())
])
pipe.fit(X_train, y_train)
# 評価
y_pred_train = pipe.predict(X_train)
y_pred_test = pipe.predict(X_test)
print(f"訓練 R²: {r2_score(y_train, y_pred_train):.4f}")
print(f"テスト R²: {r2_score(y_test, y_pred_test):.4f}")
print(f"テスト RMSE: {np.sqrt(mean_squared_error(y_test, y_pred_test)):.4f}")
# 交差検証
cv_scores = cross_val_score(pipe, X, y, cv=5, scoring="r2")
print(f"\n5-fold CV R²: {cv_scores.mean():.4f} (±{cv_scores.std():.4f})")
また,今回は目的変数が連続数なので,予測値と実測値を比較するプロットを作成してみます.全てがy=xの赤い点線上に乗れば理想です.
# 予測 vs 実測プロット
plt.figure(figsize=(7, 7))
plt.scatter(y_test, y_pred_test, alpha=0.5, s=20)
plt.plot([-12, 2], [-12, 2], "r--", label="ideal") #logSは-12~2の間に分布
plt.xlabel("Measured logS")
plt.ylabel("Predicted logS")
plt.legend()
plt.show()

R² = 0.676と,元論文の結果(R² = 0.811)
を下回っていますが、元論文は LogP を含む4記述子を使っているのに対し,今回はまだ LogP を含めていないことが主な要因と考えられます。
y=xの直線とある程度の相関がみられますが,まだ乗っていないプロットも多く改善の余地があると思われます.
Top5・Worst5分子
ちなみに,趣旨とは外れますが,このデータセットのうち溶けやすい/溶けにくいTop5を描画してみようと思います.
.nlargestとnsmallestを使うことで,抜き出すことができます.
iterrows() は各行を (インデックス, 行データ) で返します.
_ はインデックス(使わないので捨てる),row は1行分のデータ(辞書的にアクセスできる)になります.
# 溶けやすい Top 5
top5_soluble = df.nlargest(5, target_col)
# 溶けにくい Top 5
top5_insoluble = df.nsmallest(5, target_col)
# 溶けやすい分子を描画
print("=== 溶けやすい Top 5 ===")
for _, row in top5_soluble.iterrows():
print(f" {row['Compound ID']:20s}: logS = {row[target_col]:.2f}")
mols_sol = [Chem.MolFromSmiles(smi) for smi in top5_soluble["smiles"]]
legends_sol = [f"{row['Compound ID']}\nlogS={row[target_col]:.1f}"
for _, row in top5_soluble.iterrows()]
img_sol = Draw.MolsToGridImage(mols_sol, molsPerRow=5, subImgSize=(250, 200),
legends=legends_sol)
display(img_sol)
# 溶けにくい分子を描画
print("\n=== 溶けにくい Top 5 ===")
for _, row in top5_insoluble.iterrows():
print(f" {row['Compound ID']:20s}: logS = {row[target_col]:.2f}")
mols_insol = [Chem.MolFromSmiles(smi) for smi in top5_insoluble["smiles"]]
legends_insol = [f"{row['Compound ID']}\nlogS={row[target_col]:.1f}"
for _, row in top5_insoluble.iterrows()]
img_insol = Draw.MolsToGridImage(mols_insol, molsPerRow=5, subImgSize=(250, 200),
legends=legends_insol)
display(img_insol)

小さくて極性が高い → 溶ける、大きくて疎水的 → 溶けないといった傾向が顕著に現れており,
Top5/Worst5どちらも納得感のある分子だなと思います.
まとめ
今回は Delaney データセットの読み込みとEDA、6つの既存記述子による線形回帰のベースラインモデルを作成しました.
交差検証の結果はR²=0.676 と,溶解度の変動の約 68% を説明できましたが,まだ改善の余地があります.
次回は RDKit を使って記述子を自分で計算・追加し,精度の改善を狙います.




