随机划分会让分子机器学习的指标系统性虚高——同一系列的类似物被同时分到训练集和测试集,模型只要记住这个系列就能拿高分。而真实项目面对的是新骨架,这才是模型能力的真实考验。骨架划分就是为此设计的。
Bemis-Murcko 骨架
from rdkit import Chem
from rdkit.Chem.Scaffolds import MurckoScaffold
def get_scaffold(smiles, include_chirality=False):
"""提取 Bemis-Murcko 骨架:去掉侧链,保留环系与连接它们的链"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
return MurckoScaffold.MurckoScaffoldSmiles(
mol=mol, includeChirality=include_chirality)
# 例子
for s in ["CC(=O)Nc1ccc(O)cc1", # 对乙酰氨基酚
"CCC(=O)Nc1ccc(OC)cc1"]: # 类似物
print(s, "→", get_scaffold(s))
# 两者骨架相同(c1ccccc1),会被分到同一组
实现骨架划分
from collections import defaultdict
import numpy as np
def scaffold_split(smiles_list, frac_train=0.8, frac_valid=0.1,
balanced=True, seed=42):
"""balanced=True: 大骨架组优先进训练集(DeepChem 风格)
balanced=False: 严格按骨架组大小降序分配"""
scaffolds = defaultdict(list)
for i, smi in enumerate(smiles_list):
sc = get_scaffold(smi)
scaffolds[sc if sc is not None else f"__invalid_{i}"].append(i)
groups = list(scaffolds.values())
if balanced:
rng = np.random.RandomState(seed)
big = [g for g in groups if len(g) > len(smiles_list) * frac_valid / 2]
small = [g for g in groups if g not in big]
rng.shuffle(big)
rng.shuffle(small)
groups = big + small # 大组先进训练集
else:
groups = sorted(groups, key=len, reverse=True)
n = len(smiles_list)
n_train, n_valid = int(frac_train * n), int(frac_valid * n)
train, valid, test = [], [], []
for g in groups:
if len(train) + len(g) <= n_train:
train += g
elif len(valid) + len(g) <= n_valid:
valid += g
else:
test += g
return train, valid, test
tr, va, te = scaffold_split(df["clean_smiles"].tolist())
print(f"训练 {len(tr)} / 验证 {len(va)} / 测试 {len(te)}")
# 验证:训练集与测试集的骨架不应有交集
tr_sc = {get_scaffold(df['clean_smiles'].iloc[i]) for i in tr}
te_sc = {get_scaffold(df['clean_smiles'].iloc[i]) for i in te}
print("骨架交集:", len(tr_sc & te_sc), "(应为 0)")
亲手做一次对照实验
这是最值得花二十分钟做的实验——亲眼看到指标掉下来,比读十遍警告都管用:
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, r2_score
import numpy as np
def evaluate(split_fn, name, n_seeds=5):
rmses, r2s = [], []
for seed in range(n_seeds):
tr, va, te = split_fn(seed)
model = RandomForestRegressor(n_estimators=500, n_jobs=-1,
random_state=seed)
model.fit(X[tr], y[tr])
pred = model.predict(X[te])
rmses.append(np.sqrt(mean_squared_error(y[te], pred)))
r2s.append(r2_score(y[te], pred))
print(f"{name:12s} RMSE = {np.mean(rmses):.3f} ± {np.std(rmses):.3f} "
f"R² = {np.mean(r2s):.3f} ± {np.std(r2s):.3f}")
def random_split(seed):
rng = np.random.RandomState(seed)
idx = rng.permutation(len(y))
n_tr, n_va = int(0.8*len(y)), int(0.1*len(y))
return idx[:n_tr], idx[n_tr:n_tr+n_va], idx[n_tr+n_va:]
evaluate(random_split, "随机划分")
evaluate(lambda s: scaffold_split(smiles, seed=s), "骨架划分")
# 典型结果:
# 随机划分 RMSE = 0.62 ± 0.03 R² = 0.71 ± 0.02
# 骨架划分 RMSE = 0.89 ± 0.08 R² = 0.38 ± 0.09
# 差距就是「记住系列」带来的虚高
更严格的划分方式
| 划分 | 严格度 | 对应场景 |
|---|---|---|
| 随机 | 最松 | 只作为性能上界参考 |
| 骨架(Murcko) | 中 | 默认选择,遇到新骨架 |
| 通用骨架(去原子类型) | 较严 | 更彻底的骨架分离 |
| 聚类划分(Butina) | 较严 | 按整体相似性分离 |
| 时间划分 | 最严 | 最接近真实前瞻预测 |
# 通用骨架:进一步去掉原子类型信息
from rdkit.Chem.Scaffolds.MurckoScaffold import MakeScaffoldGeneric
def generic_scaffold(smi):
mol = Chem.MolFromSmiles(smi)
if mol is None:
return None
sc = MurckoScaffold.GetScaffoldForMol(mol)
return Chem.MolToSmiles(MakeScaffoldGeneric(sc))
# 时间划分:如果数据带文献年份,这是最有说服力的
df_sorted = df.sort_values("year")
n = len(df_sorted)
train = df_sorted.iloc[:int(0.8*n)]
test = df_sorted.iloc[int(0.8*n):]
# 用早期数据预测后期化合物 —— 完全模拟真实处境
报告结果的规范
- 必须写明划分方式。只报一个 RMSE 而不说划分方式,等于没报。
- 必须多种子:骨架划分本身有随机性(组的分配顺序),跑 5 个种子报均值 ± 标准差。
- 建议同时报随机划分:两者的差距本身是关于模型泛化能力的重要信息。
- 小数据集要特别谨慎:几百个分子时,不同种子间波动可能大于模型间差异。
常见坑与提示
- 随机划分的指标系统性虚高,只能作为上界参考;
- 亲手跑一次随机 vs 骨架的对照实验,胜过任何说教;
- 骨架划分也有随机性,必须多种子报均值 ± 标准差;
- 有时间信息时,时间划分最接近真实前瞻场景。