骨架划分是分子机器学习中最重要的评测实践之一。它解决的问题很具体:药物数据是成系列产生的,随机划分会让测试集分子的「近亲」留在训练集中——模型只需查表插值就能得高分。
问题的根源
# 【药物数据的产生方式】:
# 化学家围绕一个骨架合成几十到几百个类似物
# → 这些分子彼此高度相似
# → 【数据点不是独立同分布的】
#
# 【随机划分的后果】:
# 同一系列的分子被随机分到训练集与测试集
# → 测试集中的每个分子,训练集中都有它的近亲
# → 模型只需「找到最相似的训练分子,输出它的活性」
# → 【这不是泛化,是记忆】
#
# 【典型的性能差距】:
# 同一模型、同一数据:
# 随机划分 R² = 0.85
# 骨架划分 R² = 0.55
# 时间划分 R² = 0.40
#
# 【而实际部署时面对的是什么】:
# 新合成的分子(可能是新系列)
# → 更接近骨架划分或时间划分的情形
# → 【随机划分的数字对实际决策毫无参考价值】
#
# 【这个问题有多普遍】:
# 相当比例的已发表工作仍以随机划分为主要结果
# → 【读论文时必须查这一点】(见 169)
实现
from rdkit import Chem
from rdkit.Chem.Scaffolds import MurckoScaffold
from collections import defaultdict
import numpy as np
def scaffold_split(smiles_list, frac_train=0.8, frac_valid=0.1,
generic=False, include_chirality=False, seed=0):
"""按 Bemis-Murcko 骨架划分数据集"""
groups = defaultdict(list)
for i, smi in enumerate(smiles_list):
mol = Chem.MolFromSmiles(smi)
if mol is None:
continue
scaf = MurckoScaffold.GetScaffoldForMol(mol)
if generic:
try:
scaf = MurckoScaffold.MakeScaffoldGeneric(scaf)
except Exception:
pass
key = Chem.MolToSmiles(scaf)
# 【无环分子单独处理,避免全部归为一类】
if not key:
key = f"__acyclic_{i}__"
groups[key].append(i)
# 【大骨架优先进训练集】
# → 让测试集包含更多罕见骨架,测试更严格
sets = sorted(groups.values(), key=len, reverse=True)
n = len(smiles_list)
n_train, n_valid = frac_train * n, frac_valid * n
train, valid, test = [], [], []
for s in sets:
if len(train) + len(s) <= n_train:
train += s
elif len(valid) + len(s) <= n_valid:
valid += s
else:
test += s
return train, valid, test
# 【变体:随机骨架划分】
# 把骨架组随机分配,而非按大小排序
# → 可以跑多个种子,得到方差估计
def random_scaffold_split(smiles_list, frac_train=0.8, seed=0):
groups = defaultdict(list)
for i, smi in enumerate(smiles_list):
mol = Chem.MolFromSmiles(smi)
if mol is None:
continue
key = MurckoScaffold.MurckoScaffoldSmiles(mol=mol) or f"__acyclic_{i}__"
groups[key].append(i)
rng = np.random.default_rng(seed)
sets = list(groups.values())
rng.shuffle(sets)
n_train = frac_train * len(smiles_list)
train, test = [], []
for s in sets:
(train if len(train) + len(s) <= n_train else test).extend(s)
return train, test
# 【推荐用这个】:能跑多个种子,报告均值 ± 标准差
验证划分是否足够严格
# 【骨架划分不一定就够严格】
# 不同的骨架可能仍然很相似
# (如苯并咪唑 vs 苯并噻唑)
#
# 【必做的检查】:看测试集与训练集的相似度分布
from rdkit import DataStructs
from rdkit.Chem import rdFingerprintGenerator
gen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048)
def split_severity(train_smiles, test_smiles):
tr_fps = [gen.GetFingerprint(Chem.MolFromSmiles(s))
for s in train_smiles if Chem.MolFromSmiles(s)]
max_sims = []
for s in test_smiles:
m = Chem.MolFromSmiles(s)
if m is None:
continue
sims = DataStructs.BulkTanimotoSimilarity(gen.GetFingerprint(m), tr_fps)
max_sims.append(max(sims))
max_sims = np.array(max_sims)
print("测试集分子与训练集的最大相似度:")
print(f" 中位数 {np.median(max_sims):.3f}")
print(f" > 0.7 比例 {(max_sims > 0.7).mean():.1%} ← 【应该很低】")
print(f" > 0.5 比例 {(max_sims > 0.5).mean():.1%}")
print(f" < 0.3 比例 {(max_sims < 0.3).mean():.1%} ← 真正的外推")
return max_sims
# 【判读】:
# > 0.7 的比例超过 20% → 【划分不够严格】
# → 考虑用 generic 骨架或聚类划分
#
# 【更有信息量的报告方式】:
# 按相似度区间分层报告性能
def stratified_performance(y_true, y_pred, max_sims):
bins = [(0.0, 0.3), (0.3, 0.5), (0.5, 0.7), (0.7, 1.0)]
print("\n按与训练集的相似度分层:")
for lo, hi in bins:
mask = (max_sims >= lo) & (max_sims < hi)
if mask.sum() < 10:
continue
rmse = np.sqrt(((y_pred[mask] - y_true[mask]) ** 2).mean())
print(f" 相似度 [{lo}, {hi}): n={mask.sum():4d}, RMSE={rmse:.3f}")
# 【这张表直接回答】:
# 「模型在距离训练集多远的地方还可信?」
# → 【这是最有实用价值的评测输出】
各种划分方式的严格程度
| 方式 | 严格度 | 对应场景 | 需要 |
|---|---|---|---|
| 随机 | 最松 | 几乎不对应任何真实场景 | — |
| Murcko 骨架 | 中 | 遇到新化学系列 | — |
| Generic 骨架 | 较严 | 杂环变体也算同类 | — |
| 指纹聚类 | 严 | 整个相似子空间未见过 | 见 050《分子相似性聚类》 |
| 时间划分 | 最接近真实 | 用历史预测未来 | 时间戳 |
| 跨项目/跨靶点 | 最严 | 模型迁移 | 多项目数据 |
时间划分:最诚实的评测
# 【如果数据有时间戳,就用时间划分】
def temporal_split(df, date_col="date", cutoff="2024-01-01"):
train = df[df[date_col] < cutoff]
test = df[df[date_col] >= cutoff]
print(f"训练 {len(train)} (至 {cutoff}),测试 {len(test)}")
return train, test
# 【为什么最诚实】:
# 它完全模拟了生产环境:
# 用截至今天的数据训练,
# 预测明天合成的分子
# → 【自然地包含了化学空间的漂移】
#
# 【滚动评测:更完整的图景】
def rolling_evaluation(df, date_col="date", n_folds=5):
df = df.sort_values(date_col)
fold_size = len(df) // (n_folds + 1)
results = []
for k in range(1, n_folds + 1):
train = df.iloc[:k * fold_size]
test = df.iloc[k * fold_size:(k + 1) * fold_size]
# ... 训练与评测 ...
results.append({"fold": k, "n_train": len(train),
"n_test": len(test)})
return results
# 【能看到「数据越多模型越好吗」这个关键问题的答案】
# 【如果内部数据没有时间戳】:
# → 【现在就开始记录】
# 合成日期、测定日期、录入日期
# → 这是低成本高回报的数据治理
# → 【一年后你会感谢自己】
报告评测结果的正确方式
# 【一份完整的评测报告应该包含】:
#
# 1) 【多种划分的结果】
# 随机划分(作为上界参照)
# 骨架划分(主要结果)
# 时间划分(如果有时间戳)
#
# 2) 【多个随机种子】
# 均值 ± 标准差
# → 【提升必须超过 2 倍标准差才算显著】
#
# 3) 【测试集与训练集的相似度分布】
#
# 4) 【分相似度区间的性能】
# ← 【最有信息量的部分】
#
# 5) 【与强基线的比较】
# ECFP + LightGBM,同等调优预算(见 131)
#
# 6) 【适用域的界定】
# 模型在什么范围内可信
#
# 【一个反面例子】:
# 「我们的模型在 XX 数据集上达到 R² = 0.92」
# → 缺少:什么划分?几个种子?基线是多少?
# → 【这个数字无法解读】
#
# 【一个正面例子】:
# 「骨架划分(5 个随机种子),R² = 0.55 ± 0.04;
# ECFP+LightGBM 基线 0.52 ± 0.05;
# 在与训练集相似度 < 0.4 的测试分子上,
# R² 降至 0.31。
# 因此模型适用于与已有系列相似度 > 0.4 的分子。」
# → 【这个结论可以直接支撑决策】
关键要点
- 药物数据是成系列产生的,不是独立同分布的——这是随机划分失效的根源;
- 无环分子的骨架为空,必须单独处理否则会全部落在同一边;
- 骨架划分后要检查测试集与训练集的最大相似度分布,确认是否足够严格;
- 分相似度区间报告性能是最有实用价值的输出——它界定了模型的适用域。
延伸资源
- 骨架:040《Scaffold 骨架》、041《Bemis–Murcko Scaffold》;聚类划分:050《分子相似性聚类》;
- 模型外推性:167《模型外推性》;Benchmark 陷阱:169《Benchmark 陷阱》;QSAR:044《QSAR 入门》。