はじめに

前回から，有機分子の溶解度予測に取り組んでいます．2004年のESOL論文を参考に，前回はデータセットの読み込み，データ内容の確認，6つの与えられた記述子からの単純な線形回帰を行い，R2 = 0.676という結果でした．
今回はRDKit記述子を自分で追加して，同じ線形回帰で再評価したいと思います．

前回使った記述子

Minimum Degree: 最小結合次数
Molecular Weight: 分子量
Number of H-Bond Donors: 水素結合ドナー数
Number of Rings: 環の数
Number of Rotatable Bonds: 回転可能結合数
Polar Surface Area: 極性表面積

RDKit で記述子を追加する

前回作成したデータフレームを元に，まず，SMILES から Mol オブジェクトを作成してmol列を作ります．
ちゃんと全部変換できたか，mol列に対して.isnull().sum()を行い確認します．

df["mol"] = df["smiles"].apply(Chem.MolFromSmiles)
print(f"変換失敗: {df['mol'].isnull().sum()}件")

次に，RDKit で記述子を計算します．
相関がありそうなものを3つまずは加えてみます．

MolLogP：分配係数．分子が水と油(オクタノール)のどちらに溶けやすいかを表す値．LogP が高い → 油に溶けやすい（疎水的）→ 水に溶けにくい．
NumHAcceptors：水素結合アクセプター数．
NumAromaticRings：芳香環の数．既存の記述子Number of Ringsは芳香環と非芳香環を区別していなかった．

df["MolLogP"] = df["mol"].apply(Descriptors.MolLogP)
df["NumHAcceptors"] = df["mol"].apply(Descriptors.NumHAcceptors)
df["NumAromaticRings"] = df["mol"].apply(Descriptors.NumAromaticRings)

既存6 + RDKit 3 = 9記述子を使って，新しいデータフレーム（X_v2）を作ります．

desc_cols_v2 = desc_cols + ["MolLogP", "NumHAcceptors", "NumAromaticRings"]
X_v2 = df[desc_cols_v2]

モデル作成・評価

# 訓練/テスト分割
X_train2, X_test2, y_train2, y_test2 = train_test_split(
    X_v2, y, test_size=0.2, random_state=42
)

# パイプライン(標準化 + 線形回帰)
pipe2 = Pipeline([
    ("scaler", StandardScaler()),
    ("model", LinearRegression())
])
pipe2.fit(X_train2, y_train2)

# 評価
y_pred_test2 = pipe2.predict(X_test2)

print(f"テスト R²: {r2_score(y_test2, y_pred_test2):.4f}")
print(f"テスト RMSE: {np.sqrt(mean_squared_error(y_test2, y_pred_test2)):.4f}")

cross validationを行い，どれだけ改善されたか確認します．

# 交差検証
cv_scores2 = cross_val_score(pipe2, X_v2, y, cv=5, scoring="r2")

# 改善量
print(f"ステップ1: CV R² = 0.6760")
print(f"ステップ2: CV R² = {cv_scores2.mean():.4f}")
print(f"改善: {cv_scores2.mean() - 0.6760:+.4f}")

0.7832と前回のR2=0.6760から大きく改善されました．

結果の考察・まとめ

改善の理由を探るべく，各記述子と目的変数の相関を確認します．

for col in ["MolLogP", "NumHAcceptors", "NumAromaticRings"]:
    corr = df[col].corr(df[target_col])
    print(f"  {col:22s}: r = {corr:+.3f}")

MolLogPは-0.828と，全ての記述子の中で一番相関が強いことがわかります．
また，NumAromaticRingsもそこそこの相関でした．
ドメイン知識に基づいて，適切な記述子を追加することが非常に効果的であったとわかりました．
（水への溶けやすさLogSを計算するのに油への溶けやすさLogPを使うのは少しズルい気もするし，正直面白くないなとは思います）
RDKitにはたくさんの記述子があり（200個以上），今回は AI に提案してもらった3つを採用しましたが，選定の根拠は「溶解度に直接関わる疎水性の指標（LogP）が欠けている」「既存の記述子を補完する情報（アクセプター数，芳香環数）」という化学的な判断に基づいています．
実際の研究においても，知識と経験からいかに適切なものを選ぶのかが重要だと感じます．

次回は、モデルを線形回帰から非線形モデル（LightGBM）に切り替えて，さらなる改善を狙います．

他のRDKit記述子

今回は試していませんが，他に水溶解度を向上させる上で有効かもしれない記述子を10個紹介します．RDKit記述子について幅広く知っておくことは他の学習モデルを作成するときの参考になるかもしれないので，将来の自分のためのメモでもあります．

MolMR(分子屈折率)
分子の大きさと分極性を反映。MolWt と似ているが電子的性質も含む。
FractionCSP3(sp3 炭素の割合)
立体的な複雑さの指標。sp3 が多い分子は結晶パッキングしにくく溶けやすい傾向。
HeavyAtomCount(重原子数)
水素以外の原子数。分子の骨格の大きさを直接反映。
LabuteASA(Labute の近似溶媒接触表面積)
溶媒と接触する分子表面の面積。溶解の物理過程に直接関わる。
NumHeteroatoms(ヘテロ原子数)
N, O, S など炭素・水素以外の原子の数。極性に寄与する。
NumAliphaticRings(脂肪族環の数)
芳香環でない環の数。芳香環との区別で溶解性への影響が異なる。
NumSaturatedRings(飽和環の数)
完全に飽和した環の数。シクロヘキサンのような構造。
飽和環は芳香環より親水的。
MaxPartialCharge(最大部分電荷)
分子内で最も正に偏った原子の電荷。局所的な極性の指標。
BalabanJ(Balaban の J 指数)
分子のトポロジー(つながり方)を数値化。分岐度や環状性を反映。
BertzCT(Bertz の複雑度)
分子の構造的複雑さの指標。原子の種類と結合パターンから計算。
複雑な分子ほど溶けにくい傾向がある。