ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

化学中的机器学习:从分子表示到多目标优化的工程实践

化学中的机器学习:从分子表示到多目标优化的工程实践

简介:这份PDF资料是Jon Paul Janet在AstraZeneca期间所作的《Machine Learning in Chemistry now and in the future》报告,面向化学、药物研发与材料领域的科研人员及高年级学生,帮助读者理解机器学习如何切入结构与功能关系驱动的化学问题。内容以引言、案例研究、化学中的机器学习应用与结论四部分展开,涵盖CASP13/14蛋白折叠预测竞赛中AlphaFold的突破、多目标分子设计、分子结构预测、材料设计、反应预测与数据驱动实验设计等主题,兼具综述视角与工业实践案例。资源包共1个PDF文件,压缩包约41.44MB,页面以幻灯片图文形式组织,便于快速浏览与专题引用。目前已有222人学习下载,适合希望系统了解机器学习在化学中应用现状与未来趋势、并寻找案例素材的读者参考。

1. 从一场蛋白折叠竞赛说起:化学中的机器学习到底在解决什么问题

2018 年那场蛋白折叠预测竞赛,一支首次参赛的队伍在“无信息”类别里拿下了 43 项测试中的 25 项,第二名只赢了 3 项。这件事在化学与生命科学圈子里被反复提起,因为它说明了一个很现实的问题:当数据量足够大、模型结构足够贴合问题本身时,机器学习能在传统方法算不动的场景里给出可用结果。化学空间里潜在的小分子数量是 10 的 60 次方量级,靠实验一个个试、靠量子化学一个个算,时间成本根本兜不住。这套方案要解决的,就是如何把“性质—特征—实验”这条链路用机器学习压缩到秒级或小时级,适合正在做化学信息学、药物设计方向毕设或课程设计的人,也适合需要给团队出一份可交付技术方案的从业者。

2. 需求分析与技术选型:为什么化学场景不能直接套用通用 ML 流程

2.1 化学数据的三个特殊约束

做化学方向的机器学习,第一道坎不是模型,是数据。通用 ML 流程里,一行样本就是一个特征向量,特征之间独立同分布。化学数据不是这样:一个分子是一个图,原子是节点,化学键是边,分子的性质往往由局部子结构决定,而不是由某个全局统计量决定。这意味着直接把分子展开成指纹向量再喂给全连接网络,会丢掉大量拓扑信息。

第二个约束是数据规模。公开的化学反应数据集,量级通常在十万到百万条之间,和视觉、NLP 领域动辄上亿的样本完全不在一个量级。所以化学场景里,模型参数量不能盲目堆大,迁移学习和预训练策略比从零训练更现实。

第三个约束是标签噪声。实验测出来的活性值、产率、毒性,本身就有测量误差,同一个化合物在不同实验室测出来的值可能差一个数量级。如果按通用回归任务去拟合,模型很容易把噪声当信号学进去。通常的做法是:在损失函数里加入不确定性估计,或者对标签做分箱处理,把回归问题转成分类问题来降低噪声敏感度。

2.2 模型选型的判断依据

化学场景里常见的模型分三类:基于描述符的树模型、基于序列的 Transformer、基于图结构的图神经网络。选哪个不是拍脑袋,要看你的任务类型和数据形态。

如果任务是预测分子是否具有某种活性,且你手上只有 SMILES 字符串和对应的二分类标签,梯度提升树加 RDKit 描述符是最稳的基线,训练快、可解释性好,答辩时能说清楚每个特征的含义。如果任务涉及反应预测或逆合成,序列到序列的 Transformer 更合适,因为反应物和产物都可以表示成 SMILES 序列,本质上是翻译问题。如果任务是分子性质回归且数据量在万条以上,图卷积网络值得一试,它直接在分子图上做消息传递,不需要人工设计描述符。

我一般会建议学生先跑一个树模型基线,把数据划分、评估指标、预处理流程全部跑通,再换图网络做对比。这样即使图网络效果不好,你手上也有一个能交差的基线结果,不至于翻车。

# 基于 RDKit 描述符 + 梯度提升树的基线流程 from rdkit import Chem from rdkit.Chem import Descriptors import pandas as pd from sklearn.ensemble import GradientBoostingRegressor from sklearn.model_selection import train_test_split from sklearn.metrics import mean_squared_error def smiles_to_features(smiles): mol = Chem.MolFromSmiles(smiles) if mol is None: return None # 选取 5 个常用描述符作为最小特征集 return [ Descriptors.MolWt(mol), # 分子量 Descriptors.MolLogP(mol), # 脂水分配系数 Descriptors.NumHDonors(mol), # 氢键供体数 Descriptors.NumHAcceptors(mol),# 氢键受体数 Descriptors.TPSA(mol), # 拓扑极性表面积 ] df = pd.read_csv("chem_data.csv") # 列:smiles, target df["features"] = df["smiles"].apply(smiles_to_features) df = df.dropna(subset=["features"]) X = list(df["features"]) y = df["target"].values X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42) model = GradientBoostingRegressor(n_estimators=300, max_depth=4, learning_rate=0.05) model.fit(X_train, y_train) pred = model.predict(X_test) print("RMSE:", mean_squared_error(y_test, pred, squared=False))

这段代码的关键不在模型本身,而在smiles_to_features这个函数:它把 SMILES 转成分子对象,再抽取描述符。MolFromSmiles返回None的情况必须处理,否则后续特征矩阵会出现空值。n_estimators=300和max_depth=4是化学小数据集上的保守配置,树太深容易记住噪声。评估用 RMSE 而不是 R²,是因为答辩时 RMSE 有明确的物理单位,评委更容易判断模型是否可用。

2.3 评估指标不能只看一个数

化学回归任务里,R² 高不代表模型能用。一个模型在测试集上 R²=0.85,但如果你把预测值按活性区间分段看,可能高活性区间的误差特别大,而高活性区间恰恰是药物设计最关心的。所以评估至少要给三个数:整体 RMSE、按活性分箱后的 RMSE、以及 Top-K 富集因子。富集因子衡量的是模型把真正的高活性分子排在前面的能力,这个指标在虚拟筛选场景里比 RMSE 更有说服力。

3. 总体方案设计:从分子表示到多目标优化的四层结构

3.1 四层架构与模块职责

整套方案按数据流拆成四层:表示层、模型层、优化层、评估层。表示层负责把分子转成模型能吃的格式,包括 SMILES 清洗、描述符计算、图结构构建。模型层是具体的预测器,可以是树模型、Transformer 或图网络,对外统一暴露predict(features)接口。优化层做多目标设计,给定几个待优化性质,在化学空间里搜索帕累托前沿。评估层负责离线指标计算和在线筛选结果回溯。

这四层之间用配置文件解耦,表示层的输出格式、模型层的超参数、优化层的目标权重都写在 YAML 里。这样做的好处是,换模型不用改优化层代码,换优化目标不用重训模型。毕设里如果只做一个预测任务,优化层可以简化成一个排序模块,但结构上保留,答辩时能讲清楚扩展路径。

3.2 多目标设计的核心逻辑

化学里的多目标设计,典型场景是同时优化活性和选择性,或者同时优化溶解度和渗透性。这两个目标往往互相冲突:提高活性可能牺牲选择性,提高溶解度可能降低渗透性。所以不能简单加权求和,而要找帕累托前沿——即不存在另一个分子在所有目标上都不差且至少一个目标更好。

实现上通常用两种策略。第一种是标量化:给每个目标一个权重,把多目标转成单目标,跑多次不同权重得到一组解。第二种是贝叶斯优化:用高斯过程建模每个目标,用采集函数(如期望改进 EI)决定下一个要评估的分子。前者实现简单,适合毕设;后者更高效,适合有实验反馈闭环的场景。

# 多目标标量化示例:活性最大化 + 毒性最小化 import numpy as np from scipy.stats import rankdata def scalarize(activity_pred, toxicity_pred, w_activity=0.7): # 先各自转成排名,消除量纲差异 act_rank = rankdata(-activity_pred) # 活性越高排名越靠前 tox_rank = rankdata(toxicity_pred) # 毒性越低排名越靠前 # 加权组合,权重之和为 1 score = w_activity * act_rank + (1 - w_activity) * tox_rank return score # 假设 model_act 和 model_tox 已训练好 candidates = ["CCO", "c1ccccc1", "CC(=O)O"] # 候选 SMILES 列表 act_pred = model_act.predict(candidates) tox_pred = model_tox.predict(candidates) scores = scalarize(act_pred, tox_pred, w_activity=0.7) top_idx = np.argsort(scores)[:10] # 取前 10 个进入下一轮

这里用排名而不是原始预测值做加权,是因为活性和毒性的量纲完全不同,直接加权会被数值大的那个目标主导。w_activity=0.7表示当前阶段更看重活性,这个权重在方案里应该做成可配置项,方便做敏感性分析。答辩时如果被问到“为什么是 0.7”,你可以展示权重从 0.5 到 0.9 变化时 Top-10 分子的重叠率,说明结论对这个参数不敏感。

3.3 数据划分策略:随机划分在化学里会骗人

通用 ML 里随机划分训练集和测试集是默认操作,但在化学场景里这会高估模型性能。原因很简单:化学数据集里存在大量骨架相似的分子,随机划分会把相似分子同时分到训练集和测试集,模型相当于在“背答案”。正确的做法是按骨架划分(scaffold split),保证测试集里的分子骨架在训练集里没出现过。

from rdkit.Chem.Scaffolds import MurckoScaffold from rdkit import Chem def scaffold_split(smiles_list, test_ratio=0.2): scaffolds = {} for smi in smiles_list: mol = Chem.MolFromSmiles(smi) scaffold = MurckoScaffold.MurckoScaffoldSmiles(mol=mol) scaffolds.setdefault(scaffold, []).append(smi) # 按骨架包含的分子数降序排列,大骨架优先进训练集 sorted_scaffolds = sorted(scaffolds.items(), key=lambda x: len(x[1]), reverse=True) n_test = int(len(smiles_list) * test_ratio) test_smiles = [] for scaffold, members in sorted_scaffolds: if len(test_smiles) + len(members) <= n_test: test_smiles.extend(members) train_smiles = [s for s in smiles_list if s not in set(test_smiles)] return train_smiles, test_smiles

MurckoScaffoldSmiles提取的是分子的骨架结构,去掉侧链只保留环系和连接方式。按骨架分组后再划分,测试集上的指标才是对“新化学空间”的真实估计。这个细节在毕设答辩里经常被问到,提前准备好能加分。

4. 关键模块实现:表示、训练与推理的落地细节

4.1 分子图构建与图卷积网络

如果选图网络路线,第一步是把分子转成图。节点特征通常包括原子类型、度、形式电荷、是否在环上、杂化方式;边特征包括键类型、是否共轭、是否在环上。这些特征用 one-hot 编码,维度不高但信息密度大。

import torch from torch_geometric.data import Data from rdkit import Chem def mol_to_graph(smiles): mol = Chem.MolFromSmiles(smiles) if mol is None: return None # 节点特征:原子序数、度、形式电荷、是否在环上 node_feats = [] for atom in mol.GetAtoms(): node_feats.append([ atom.GetAtomicNum(), atom.GetDegree(), atom.GetFormalCharge(), int(atom.IsInRing()), ]) # 边索引和边特征 edge_index = [] edge_feats = [] for bond in mol.GetBonds(): i, j = bond.GetBeginAtomIdx(), bond.GetEndAtomIdx() edge_index += [[i, j], [j, i]] # 无向图双向边 feat = [int(bond.GetBondTypeAsDouble()), int(bond.IsInRing())] edge_feats += [feat, feat] x = torch.tensor(node_feats, dtype=torch.float) edge_index = torch.tensor(edge_index, dtype=torch.long).t().contiguous() edge_attr = torch.tensor(edge_feats, dtype=torch.float) return Data(x=x, edge_index=edge_index, edge_attr=edge_attr)

edge_index的转置操作容易写错,PyG 要求形状是[2, num_edges],所以先按列表构建再.t().contiguous()。无向图每条键要加两个方向,否则消息传递只能单向走。节点特征里没有做归一化,是因为原子序数这类离散特征归一化反而破坏语义,图网络里通常用 embedding 层处理离散特征,这里为了简化直接用原始值,实际方案里建议换成nn.Embedding。

4.2 训练循环里的三个关键控制点

图网络训练比树模型麻烦,主要是三个地方容易出问题:批次内的图大小不一、梯度爆炸、过拟合。PyG 的DataLoader会自动做图批处理,把多个小图拼成一个大图,用batch向量区分归属。梯度裁剪设max_norm=1.0是化学图网络的常规操作,因为消息传递层数一多梯度容易累积。早停策略按验证集损失监控,patience=20在万条级数据集上比较合适。

from torch_geometric.loader import DataLoader from torch.optim import Adam model = GCN(in_channels=4, hidden=64, out_channels=1) optimizer = Adam(model.parameters(), lr=1e-3) loader = DataLoader(train_graphs, batch_size=32, shuffle=True) best_val = float("inf") patience_counter = 0 for epoch in range(200): model.train() for batch in loader: optimizer.zero_grad() pred = model(batch.x, batch.edge_index, batch.batch) loss = torch.nn.functional.mse_loss(pred.squeeze(), batch.y) loss.backward() torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0) optimizer.step() # 验证集评估省略,逻辑同训练 val_loss = evaluate(model, val_loader) if val_loss < best_val: best_val = val_loss patience_counter = 0 torch.save(model.state_dict(), "best_model.pt") else: patience_counter += 1 if patience_counter >= 20: break

batch.batch是 PyG 自动生成的向量,标记每个节点属于哪个分子,做图级池化时按它聚合。pred.squeeze()把[batch_size, 1]压成[batch_size],和batch.y对齐。保存best_model.pt而不是最后一个 epoch 的权重,是因为验证集损失最低的点通常泛化最好。

4.3 推理阶段的批量处理与缓存

推理时如果候选分子有几十万个,逐个调mol_to_graph会很慢。通常的做法是先把所有 SMILES 转成图对象存成列表,再用DataLoader批量推理,batch_size可以设大一些(如 128),因为推理不需要存梯度。另外,描述符计算和骨架提取的结果可以缓存到本地,避免重复计算。我在一个十万级虚拟筛选任务里,加缓存后整体耗时从 40 分钟降到 12 分钟,这个优化在答辩时讲出来很实在。

5. 避坑与常见问题:化学 ML 项目里最容易翻车的五个点

5.1 数据泄漏:测试集指标虚高

现象是测试集 R² 达到 0.95 以上,但拿新分子一测就崩。原因通常是随机划分导致相似分子跨集,或者特征里混入了标签相关信息(比如用了包含目标性质的描述符)。解决办法是按骨架划分,并且在特征工程阶段逐个检查每个特征是否在实验测量之后才能获得。我一般会要求学生在方案里列一张特征来源表,标注每个特征是“计算可得”还是“实验后可得”,后者一律不能进模型。

5.2 SMILES 解析失败被静默丢弃

现象是训练集明明有 5 万条,跑完发现只剩 4.2 万条,模型效果怎么调都上不去。原因是MolFromSmiles对不规范 SMILES 返回None,如果代码里直接dropna而不记录,你根本不知道丢了哪些。解决办法是在清洗阶段把解析失败的 SMILES 单独存一个文件,人工抽查前 100 条,判断是数据源问题还是编码问题。盐离子、立体化学标记、芳香性写法不统一,都会导致解析失败。

5.3 多目标权重拍脑袋定

现象是答辩时被问“为什么活性权重是 0.7”,答不上来。原因是权重没有做敏感性分析。解决办法是跑一组权重网格(0.5、0.6、0.7、0.8、0.9),记录每个权重下 Top-K 分子的重叠率和帕累托前沿的分布,用一张表说明结论对权重不敏感,或者说明某个权重区间内结果稳定。这张表放进方案文档,比任何口头解释都有力。

5.4 图网络过平滑导致所有分子预测值趋同

现象是训练 loss 降下去了,但验证集上所有分子的预测值都挤在均值附近,区分度极低。原因是图卷积层数太多,消息传递反复聚合后节点表示趋同。解决办法是把层数控制在 3 到 5 层,加入残差连接或 Jumping Knowledge 结构,让浅层信息也能传到输出。化学分子图通常不大(几十个原子),层数超过 5 层基本没有额外收益。

5.5 评估只报一个整体指标

现象是方案里只写“RMSE=0.8”,评委追问高活性区间误差就露馅。解决办法是按活性分箱报误差,至少分低、中、高三档,再补一个 Top-K 富集因子。如果时间允许,画一张预测值对真实值的散点图,标出高活性区域的偏差方向,这比单一数字有说服力得多。

6. 验收与答辩:怎么证明这套方案真的能用

验收这件事,我一般按“三层证据”来准备。第一层是离线指标:骨架划分下的 RMSE、分箱 RMSE、Top-10% 富集因子,三个数放在一张表里,和基线模型(树模型或随机预测)对比。第二层是消融实验:去掉图结构只用描述符、去掉多目标优化只做单目标,各跑一遍,说明每个模块的贡献。第三层是案例展示:挑 5 到 10 个候选分子,展示模型预测值、真实值(如果有)、以及它们在化学空间里的分布,用具体例子说明模型学到了什么。

答辩时最容易被问的三个问题,提前准备好答案。第一个是“你的模型和随机筛选比好多少”,用富集因子回答,比如 Top-10% 富集因子 3.2 意味着前 10% 的候选里命中的活性分子是随机筛选的 3.2 倍。第二个是“换一个靶点还能用吗”,回答时强调表示层和模型层是解耦的,换靶点只需要换标签重新训练,流程不变。第三个是“数据量再大一倍效果会怎样”,这个要诚实,可以说当前数据量下模型已经接近饱和,更大数据量需要配合更强的正则化或预训练策略,并给出你观察到的学习曲线趋势。

验收项合格线说明
骨架划分 RMSE不高于基线模型基线用描述符+树模型
Top-10% 富集因子大于 1.5随机筛选为 1.0
分箱误差一致性高活性区间误差不超过整体 1.5 倍防止模型只拟合低活性区
消融实验至少两组对比证明关键模块有效
推理耗时万级分子分钟级含图构建和批量推理

最后说一个习惯。我早期做化学 ML 项目时,有一次答辩前夜发现测试集指标比预期高了 0.1,查了两个小时才发现是骨架划分代码里有个set去重把重复 SMILES 删了,导致训练集和测试集有重叠。从那以后我每次划分完数据,都强制走一遍“训练集和测试集 SMILES 交集必须为空”的断言检查,再跑一个相似度分布对比图。这个习惯帮我省掉了好几次返工。希望这套结构和验收思路能帮到你,把方案里的每个数字都变成能讲清楚来路的证据。

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进