Mendelevium
Diary
Drug Design
Field Knowledge
Academia
Yang
Biology
Physics
Free Energy
Machine Learning & AI
Active Learning
Basics
Boltz-2
Data
Generation
Interpretability
QSAR application
Representations
Mol2Image
Workflow & Agent
Molecular Dynamics
FF & Algorithm
Small Molecule
martini
water
Interaction
Modeling & Tools
QM
Sampling & Analysis
Allostery
Fundamental
Other
Specific Sytems
Enzyme Engineering
Fiber & LLPS
Membrane
orientation_penetration
Metal
Nano Polymers
Skin Permeation
Techniques
Linux
Python
Research
Web
about
Home
Contact
Copyright © 2025 Xufan Gao | Academic Research Blog
Home
>
Machine Learning & AI
> Interpretability
A Bunch of Biophysics is Loading ...
Interpretability
药物协同预测黑箱怎么破?用graphlet指纹把子结构直接送进模型
药物协同预测黑箱怎么破?用graphlet指纹把子结构直接送进模型 本文信息 标题:GraFSyn:An Interpretable Deep Learning Framework for Anticancer Drug Synergy via Graphlet Fingerprints 作者:Wei Xia, Yayu Tian, Shiyu Zhou, Huanhuan Du, Mingchen Xiao, Zhuxu Ge, Xuan He 发表期刊:Journal of Chemical Information and Modeling 发表时间:2026年(Received February 12, 2026;Revised May 9, 2026;Accepted May 13, 2026) DOI:https://doi.org/10.1021/acs.jcim.6c00458 单位:中国辽宁省沈阳市,东北大学医学院与生物信息工程学院 引用格式:Xia, W.; Tian, Y.; Zhou, S.; Du, H.; Xiao, M.; Ge, Z.; He, X. GraFSyn: An Interpretable Deep Learning Framework for Anticancer Drug Synergy via Graphlet Fingerprints. Journal of Chemical Information and Modeling. 2026. https://doi.org/10.1021/acs.jcim.6c00458 代码与数据:https://github.com/drug-XW/GraFSyn.git 摘要 预测药物协同对加速发现有效的抗癌联合疗法具有重要意义。协同效应本质上依赖于关键化学子结构在特定细胞环境中的精确相互作用。然而,当前基于分子图的计算方法通常依赖于隐式的原子级特征聚合,这可能掩盖关键化学子结构的拓扑表示,并限制结构可追踪性。因此,我们提出了基于Graphlet指纹的协同预测框架(GraFSyn),这是一个用于抗癌药物协同预测的深度学习框架,它使用graphlet指纹将药物编码为明确的连通子结构单元,保留了预定义的化学子结构及其拓扑特征。我们进一步引入了动态多尺度卷积模块(DMSC),以便从高维和稀疏的graphlet特征中学习信息丰富的表示。该框架还包括一个交互模块,用于捕捉药物子结构与细胞系基因表达之间依赖于上下文的相互作用。在Merck和AstraZeneca基准数据集上,GraFSyn分别实现了0.972/0.912和0.823/0.906的ROC-AUC/PR-AUC值,优于代表性基线方法。此外,归因信号可以映射回特定的药效团区域,支持协同相互作用的子结构级分析。总体而言,GraFSyn为抗癌药物组合筛选提供了一种准确且结构可追踪的方法。 核心结论 GraFSyn使用graphlet指纹编码药物:将药物表示为明确的连通子结构单元,保留预定义化学子结构及其拓扑特征 DMSC模块处理稀疏高维特征:从高维和稀疏的graphlet特征中学习信息丰富的表示 交互模块捕捉上下文依赖:建模药物子结构与细胞系基因表达之间的环境依赖相互作用 性能优于代表性基线方法:在Merck数据集上达到ROC-AUC/PR-AUC为0.972/0.912,在AstraZeneca数据集达到0.823/0.906 归因信号可映射到药效团区域:支持协同相互作用的子结构级分析,提供结构可追踪的药物组合筛选方法 背景 药物协同预测的结构可解释性困境 药物协同预测的核心目标:两种药合在一起,是产生加成效应、相互抵消,还是能把肿瘤细胞压得更狠?这个问题在抗癌药物研发中尤其重要,因为临床上很多有效方案都是药物联用而非单药高剂量。然而,传统深度学习方法虽然能给出不错的预测分数,却常常难以回答一个更实际的问题:模型到底看中了分子的哪一块? 当前基于分子图的 GNN 方法普遍采用隐式原子级特征聚合策略,存在几个关键问题: 结构语义在层层压缩中被磨平:通过若干层卷积、池化和读操作将原子和键的初始特征压缩成单一表示向量,端到端训练虽方便,但代价是丢失关键的子结构信息,这意味着模型难以回答“哪些子结构、哪些药效团、哪些局部化学模式在驱动协同” 难以追踪具体化学区域:对药物发现而言,真正有用的不是一串漂浮的 embedding,而是能明确回答”是因为哪个芳香环“、”因为哪类取代基“还是”因为某个特定的药效团模式“这类具体问题 解释停留在抽象层面:如果模型只能说”某个隐藏维度重要“或”某个特征向量权重更大“,这种解释对后续药物优化的帮助有限,药化人员需要的是能回到具体化学区域的指导 更关键的问题是,这类黑箱模型难以回答药物化学家和生物学家的实际问题:是因为哪个芳香环产生作用?是因为哪类取代基贡献显著?还是因为某个特定的药效团模式在驱动协同?如果模型只能说“某个隐藏维度重要”,这种解释对后续药物优化帮助有限。 Graphlet Fingerprints 的独特优势 Graphlet fingerprints 是一种显式编码预定义子结构的表示方法。与传统的 ECFP 等扩展连通性指纹不同,它不是枚举所有路径或环形子结构,而是枚举连通诱导子图并按节点数和原子/键属性进行分类计数。这意味着每个 graphlet 都有明确的化学含义:5节点的苯环、3节点的酰胺片段、包含氮和氧的特定模式,等等。 这种表示方式有几个关键优势:模型从输入层就知道自己在处理哪些子结构,而不是在黑箱里自行推断官能团;不同的 graphlet 可以组合成更复杂的模式,不像某些固定模式指纹那样僵化;药物化学家习惯从子结构和官能团角度思考问题,graphlet 正好符合这种思维模式。 当然,graphlet 也有自己的问题。最大的挑战在于稀疏性:预定义的子结构类型很多,但单个分子通常只包含其中一小部分,导致大多数 graphlet 计数为零。这种高维稀疏输入直接喂给深度模型很难训练,这也是 GraFSyn 引入 DMSC 模块的直接动机。 传统方法的可解释性局限 当前药物协同预测领域的可解释性方法大致分为三类,各有明显局限:基于注意力的方法虽然能标出原子重要性,但注意力权重往往过于分散,很难形成清晰的结构级解释;基于梯度的归因能反推哪些原子对预测贡献更大,但容易受局部梯度消失或爆炸影响,且难以捕捉全局结构模式;事后解释器(如 LIME、SHAP)在训练好的黑箱模型外再加一层解释,这种方法的问题是“预测和解释说两套话”,且特征替换未必符合化学合理性。 GraFSyn 的思路与这些都不同。它从输入定义开始就把可解释性写进架构里。graphlet fingerprints 确保了”模型看什么”和”人能读懂什么”之间的对应关系,DMSC 确保了这种结构语义能被模型有效利用。 关键科学问题 如何把药物的明确子结构信息保留下来:避免一上来就压成不可追踪的隐向量,让模型从一开始就站在官能团和局部拓扑的层面上 如何处理graphlet计数这种稀疏高维输入:原始graphlet很稀疏,直接喂给模型会不好学,需要专门设计把能看懂与能训练接起来 如何把预测结果映射回具体化学区域:让解释能落到药效团和官能团上,而不是停在抽象权重或热图上 如何建模药物与细胞背景的交互:协同效应高度依赖细胞系,模型不能只看药本身,还得看这套细胞背景下药会怎么表现 研究内容 整体架构:从 graphlet 到协同预测的完整管线 GraFSyn 的整体框架可以看作一条清晰的流水线:分子图 → graphlet fingerprints → DMSC → 特征精炼(CR) → 交互模块(IM)→ 协同预测。这条线上每个环节都有明确职责,不是简单地把组件堆在一起。 图1:GraFSyn 的整体框架。GFE 先从分子图中提取 graphlet fingerprints,DMSC 把稀疏计数变成密集表示,CR 做特征精炼,IM 负责药物子结构与细胞系基因表达的交互,最后进入协同预测。 这条管线的关键设计哲学是可解释性前置:从输入定义开始就确保模型能看到的东西是人类能理解的化学单元,而不是等到最后再拿解释器补救黑箱。 Graphlet Fingerprints:显式子结构编码 graphlet fingerprints 的构建过程可分为三个层次: 第一层是连通诱导子图枚举:给定分子图,系统性地枚举所有不超过6个节点的连通诱导子图。这里的“诱导”意味着子图继承原分子图中节点的原子属性和键属性,不是抽象拓扑图。原文中原子特征包括元素类型、形式电荷和芳香性,键特征包括单键、双键、三键和芳香键 第二层是graphlet分类:每个子图由拓扑结构、原子属性和键属性共同决定类别。1节点graphlet由原子属性定义,2节点graphlet由两个端点原子和连接键定义,更高阶graphlet通过递归哈希函数映射到唯一的同构类别 第三层是计数和归一化:统计每个graphlet同构类别在分子中出现的精确次数,形成频率直方图$f(G)$,再归一化成与分子大小弱相关的 concentration profile $z_G$。因此,模型输入的每一维都对应一类预定义、可追踪的化学子结构 图2:graphlet fingerprint 的构建流程。分子先表示成带属性的图,再枚举连通诱导子图,按节点数和原子/键属性分配类型,并把计数归一化。 与 ECFP 等传统指纹相比,graphlet fingerprints 的独特之处在于每个维度都有明确的化学含义。 ECFP 的某个哈希位可能对应多种不同子结构的混合,难以说清该位到底代表什么。graphlet 的每个计数都对应一类特定子结构,归因时可直接回到“这个芳香环”“这个杂环片段”等药化语言。 DMSC:处理稀疏高维输入的关键模块 DMSC(Dynamic Multi-Scale Convolution)是 GraFSyn 的核心创新之一。其设计动机是处理高维稀疏的 graphlet 计数输入。 在典型的药物分子数据集里,可能的 graphlet 类别数以千计,但单个分子通常只包含其中几十到几百类。如果直接把这种稀疏向量喂给全连接层,模型很难学到有效的权重表示,大多数位置在大部分样本里都是零,梯度信号很弱。 DMSC 的解决思路是多尺度卷积,它通过三种机制来处理稀疏的graphlet计数: 多尺度卷积核:使用$k = 3$、5、7的并行一维卷积核,从graphlet频率中捕捉不同粒度的局部相关和motif共现模式 动态自适应加权:DAWM模块受selective kernel networks启发,对不同卷积尺度做全局平均池化,再用softmax得到尺度权重$\alpha_k$,让模型根据具体分子在细粒度原子模式和粗粒度环系模式之间调整关注点 压缩与池化:多尺度特征加权后,还经过最大池化、$1 \times 1$压缩层和AdaptiveAvgPool1d,得到固定维度的紧凑药物表示$H_d$;细胞系基因表达特征也经过类似流程得到$H_c$ 消融实验清楚证明了 DMSC 的重要性:去掉该模块后,F1 分数和 PR-AUC 都明显下降。这表明原始 graphlet 计数虽有明确的化学含义,但直接用作特征表示过于粗糙,必须经过适当变换才能被深度模型有效利用。 交互模块:建模药物与细胞背景的依赖关系 抗癌药物协同效应的一个关键特点是高度依赖细胞背景。同一对药物在不同细胞系中的协同强度可能差异巨大,这是因为不同细胞系的基因表达谱、代谢状态、信号通路活性各不相同。 GraFSyn 的交互模块(Interaction Module)专门处理这个问题。它不是简单地把药物表示和细胞系特征拼接,而是通过学习的交互模式来捕捉药物子结构与细胞基因表达之间的依赖关系。 具体来说,IM 不是简单拼接,也不是普通点乘交互,而是使用双向 cross-attention机制: cell-guided drug view:以药物特征作为Query,细胞系特征作为Key和Value,回答“在这个细胞背景下,哪些药物子结构更关键” drug-guided cell view:反过来以细胞系特征作为Query,药物特征作为Key和Value,建模“哪些细胞表达模式可能被药物motif调制” 残差与层归一化:cross-attention输出与原表示相加后做LayerNorm,再用全局平均池化压成固定长度向量 药物对称融合:对药物A和药物B的上下文表示做逐元素加和、逐元素乘积和绝对差拼接,保证药物输入顺序不影响最终预测 该设计将协同建模从”只看药”推进到”药+细胞一起看”。许多药物协同失效或毒副作用增强,恰恰来自细胞背景差异,若模型忽略这一点,在实际应用中容易出问题。 训练和数据 GraFSyn 的 benchmark 分成两类:一类是 Merck,一类是 AstraZeneca。两者都用于抗癌药物协同预测,但规模、类别比例和化学空间不同,因此不能只看一个平均分。 Merck原始数据包含22,737个drug-cell line triplets,覆盖38个药物和39个癌细胞系。作者用Combenefit按Loewe additivity模型计算协同分数,Loewe分数大于30定义为协同样本,小于0定义为拮抗样本,0到30之间的模糊样本被排除。预处理后留下10,650个有效triplets,其中1,973个阳性、8,677个阴性,覆盖36个药物和31个细胞系 AstraZeneca用于测试更宽的化学空间。该数据集覆盖52个药物和24个癌细胞系,同样使用上述阈值后得到668个有效triplets,其中480个阳性、188个阴性。它的类别比例和Merck明显不同,所以论文在训练损失中按数据集类别分布设置权重 细胞系表达来自Cancer Cell Line Encyclopedia。作者使用LINCS L1000定义的977个landmark genes,在控制维度的同时保留广泛转录组信息。原文说这些基因约能捕捉总转录组变异的80% 药物结构来自SMILES。标准SMILES主要从PubChem获得,少数无法匹配的化合物由ChEMBL和DrugBank补充验证,并用RDKit解析与校验 训练目标是二分类协同概率,预测层为MLP加sigmoid。由于协同样本和非协同样本比例不均衡,作者使用加权二元交叉熵,按训练集类别分布重新平衡正负样本梯度贡献;文中给出的权重参数为Merck的0.818和AstraZeneca的0.285。标准比较采用一致的5-fold cross-validation,并和LR、RF、XGBoost、DeepSynergy、DTF、DeepDDS-GAT、DeepDDS-GCN、MPFFPSDC、DFFNDDS、SynergyGTN等十个基线比较。 更重要的是,论文没有只报告随机切分,还做了三类冷启动评估:Leave-Combination-Out、Leave-Cell-Line-Out和Leave-Drug-Out。LCO考察未见药物组合,LCLO考察未见细胞系,LDO考察完全未见药物;其中LDO最接近新药筛选场景,也最难。 创新点 将子结构语义直接融入输入层,避免模型在黑箱中自行推断官能团 用 DMSC 处理稀疏 graphlet 计数,将难学的离散输入转换为连续表示 将解释集成到架构设计,预测和归因使用相同的结构线索 许多“可解释模型”仅在输出端添加解释器,模型本身仍是黑箱。GraFSyn 将子结构池化、权重映射和归因回传整合到同一链条中:模型所见与最终解释尽量一致。这虽不保证解释百分之百正确,但至少减少了“预测和解释不一致”的问题。 研究结果 主要基准测试结果 在三个不同测试设置中,GraFSyn 均取得良好表现: 场景 GraFSyn 的表现 Merck ROC-AUC 0.972,PR-AUC 0.912 AstraZeneca ROC-AUC 0.823,PR-AUC 0.906 Leave-Drug-Out ROC-AUC 0.798,PR-AUC 0.526 数值本身不是最突出,真正有意义的是其在冷启动场景下的稳定表现。这对药物协同任务至关重要,因为实际应用常遇到未见新药。许多方法在训练集上表现优异但 Leave-Drug-Out 性能大幅下降,而 GraFSyn 证明将子结构语义融入输入层确实有效,模型因此保留了可迁移信息。 从数值看,Leave-Drug-Out 的 PR-AUC 0.526 虽低于主结果,但考虑到任务难度——模型预测的是包含未见药物的组合——这一下降是预期中的。更关键的是,DeepDDS-GAT在同一LDO设置下PR-AUC为0.370、F1为0.194,降幅更明显;GraFSyn仍保留了较强的未见药物泛化能力。 与基线方法的对比 GraFSyn 在三类基线方法上都有明显优势: 对比传统机器学习:如随机森林、SVM等基于手工特征的方法,通常需要专家精心设计特征,且对新结构的泛化能力有限 对比传统GNN:如GraphCNN、GAT等基于原子级卷积的方法,在标准任务上表现不错,但在需要结构可解释性的协同预测场景下,难以解释模型为何如此预测 对比基于ECFP的方法:ECFP虽也是子结构指纹,但基于路径和环的枚举,非明确的连通诱导子图,且缺乏针对稀疏性的专门设计 GraFSyn 的优势不仅体现在分数上,还体现在结果的稳定性。在 AstraZeneca 等外部数据集上,许多方法遇到分布变化时性能明显下降;GraFSyn 证明,将子结构语义融入输入层确实有效,模型因此保留了可迁移信息。这一设计思路比单纯的百分比提升更重要。 消融实验:验证各模块必要性 为证明 GraFSyn 的高分并非偶然,作者进行了三类补充验证: 特征置换控制实验:随机打乱 graphlet 特征后,性能明显下降。这直接证明模型确实在学习 graphlet 特征中的结构信息,而非仅凭细胞系表达猜测结果。若模型能通过细胞背景解决问题,药物结构信息即为多余,这显然不合理。 顺序消融实验:逐个移除 GraFSyn 模块后,DMSC 影响最大。原文报告去掉DMSC后F1下降0.046、PR-AUC下降0.036。这与主文判断一致:DMSC 非装饰性组件,而是将稀疏子结构输入转换为可学习表示的关键步骤。缺少这一步时,原始 graphlet 计数虽可解释,但过于稀疏,深度模型难以从中学习稳定表示。 图3:消融结果。作者比较完整GraFSyn、去掉DMSC、去掉CR、去掉IM,以及用2048维Morgan fingerprint替换GFE的版本。去掉DMSC后,F1和PR-AUC掉得最明显,说明原始graphlet计数虽然可解释,但太稀疏,必须靠DMSC把它接到可学习的连续表示上。 欠采样实验:药物协同数据集通常存在严重类别不平衡,协同 pair 远少于非协同 pair。简单下采样可平衡类别,但会浪费大量数据。作者测试了不同采样比例下的模型表现,发现 GraFSyn 对类别比例变化具有一定耐受性,不仅限于“理想数据分布”。 这些验证表明:GraFSyn 的高分并非偶然,其结构设计和训练流程确实有效。如果某模块仅是偶然作用,顺序消融会立即暴露;如果依赖过拟合,Leave-Drug-Out 会明显下降。GraFSyn 在这些测试中表现稳定,说明其设计思路是可靠的。 解释性结果 GraFSyn 的另一重要特点是归因可映射到具体化学区域。 图4:双向cross-attention权重热图。行对应不同drug-cell line样本,列对应不同attention方向;作者把权重聚合成$16 \times 16$特征块,横轴是药物子结构特征块,纵轴是细胞系表达特征块,颜色越深表示权重越高。 这张图回答的是模型是否真的在按上下文重排子结构重要性。原文用Cyclophosphamide在T47D中的配对案例说明:当它与Zolinza配对时,注意力分布和与BEZ-235配对时明显不同,说明GraFSyn并没有把药物表示固定死,而是会随配对药物和细胞背景改变关注区域。 原文还比较了同一药对在T47D和LNCaP中的权重分布。作者将这种差异与两类细胞已知的遗传背景差别联系起来,例如PIK3CA突变和PTEN缺失等。这一点很关键:GraFSyn 的交互模块不是只在数学上“做了attention”,而是确实学到了细胞背景改变后,哪些化学子结构更值得看。 图5:训练前后embedding的t-SNE投影。列对应CAOV3、KPL1、SW837三个细胞系,行对应训练前后;蓝点是协同样本,黄点是非协同样本。 这张图属于表示学习的 sanity check。训练前,协同和非协同样本大幅混在一起,边界很模糊;训练后,三种细胞背景下都出现了更清晰的分离趋势,只是分离程度因细胞系而异。原文据此认为,GraFSyn 确实把drug-cell interaction投影到了更有判别力的潜空间里,而不是仅仅在输出层“调阈值”。 图6:药效团热点的归因可视化。深红色表示贡献更高。5-FU/ABT-888 在不同细胞系中会把关注点从 C-F 键切换到嘧啶骨架,这说明模型学到的是上下文相关的结构信号。 作者还提供了具体案例:Vorinostat 和 Bortezomib 在 MSTO-211H 细胞系的协同预测中,模型将贡献较高的区域指向羟肟酸片段、芳香帽和蛋白酶体抑制相关结构,这与已知机制相符。这表明 GraFSyn 的解释不仅停留在热图层面,至少在典型药物对上能与已知药理逻辑对应。 该案例说明:同一对药在不同细胞系中,关注的子结构可能不同。这与真实生物环境一致。GraFSyn 输出的不是孤立“药物标签”,而是带上下文的结构解释。对药物联用研究而言,这种信息比简单高分更有价值,因为它能直接解释背景变化时模型预测改变的原因。 图7:Vorinostat-Bortezomib 的 leave-combination-out 案例。图里把单药作用、协同机制和细胞层面解释放到一起,方便读者把模型输出和已知药理机制对照起来。 该图不仅展示“模型认为重要的区域”,还将这种重要性连接到已知机制。Vorinostat 侧重 HDAC 抑制,Bortezomib 侧重蛋白酶体抑制,两者叠加指向更强的肿瘤细胞杀伤。GraFSyn 的归因结果重建了这一机制链条。 此外,这类案例图说明了:可解释性不是附加组件,而是筛选候选组合时的关键证据。如果一个模型只能预测协同分数却无法说明原因,药化人员仍需自行分析。GraFSyn 将这一步骤前移,虽非真正的机制证明,但已比纯黑箱预测更实用。 关键结论与批判性总结 优点:GraFSyn 的路径清晰,graphlet fingerprints、DMSC、交互模块和解释模块分工明确,并非简单堆砌组件 局限:graphlet 依赖预定义子结构,表达能力受子结构库约束;若真实药效来自细微构象变化,该表示可能无法完全捕捉 未来方向:将该方法应用于更多外部药物组合和不同肿瘤背景,验证其解释是否仍能与已知机制稳定对应,而非仅在 benchmark 上成立 GraFSyn 的核心价值在于推动协同预测的发展:模型能回答”为什么是这个结构””为什么在此细胞背景下是这个解释”。这类方法需要更多数据集和外部验证才能走向实用。 一句话总结:它不仅构建了更强的预测器,还将药物子结构、细胞上下文和解释结果形成完整链条。这条链条若能在更多任务上验证,才真正有价值。 补充验证 超参数设置:batch size = 32,最大图节点数 = 6,优化器 = Adam,调度器 = ReduceLROnPlateau,L2正则化 = $1 \times 10^{-5}$,预测层dropout = 0.3。该配置较为克制,非依靠极端大模型获得的结果 特征置换控制实验:打乱特征后性能下降,说明模型确实在学习graphlet特征,不是只依赖数据集偏差 顺序消融实验:移除DMSC影响最大,与主文判断一致 欠采样实验:不平衡数据下,模型仍能保持可用表现,说明其不仅限于理想数据分布 仅看主文图表可能误解为“新指纹 + GNN”。SI 提示,真正有效的是整套输入表示、稀疏处理、上下文交互和验证流程。 本文没有将可解释性停留在热图层面,而是将解释前移:从输入定义到模型池化再到案例回溯形成完整链条。对依赖上下文的药物协同任务,这种方法比单独提供 attention map 更可靠。
Machine Learning & AI
· 2026-05-31
TradePool:用PubChem指纹子结构池化与映射,给GNN分子性质预测提供可量化的原子归因
TradePool:用PubChem指纹子结构池化与映射,给GNN分子性质预测提供可量化的原子归因 本文信息 标题:TradePool:一种用于量化分子性质预测中原子归因的新型可解释框架 作者:Bingwei Ni, Wanxiang Shen(申万祥), Zhuyifan Ye* 发表时间:2025年12月22日 单位:澳门理工大学(中国澳门),宁波大学药物发现技术研究院(中国浙江),浙江大学药学院(中国杭州) 引用格式:Ni, B.; Shen, W.; Ye, Z. TradePool: A Novel Interpretable Framework for Quantifying Atomic Attribution Values in Molecular Property Prediction. J. Chem. Inf. Model. 2025, 65, XXX–XXX. https://doi.org/10.1021/acs.jcim.5c02225 开源代码与数据:https://github.com/nibingwei123/TradePool 摘要 图神经网络的可解释性一直是化合物性质预测领域的焦点。GNN在小样本化合物数据集建模上表现良好,但现有可解释方法难以准确解释原子归因值(单个原子对模型预测贡献的定量度量),使得先导化合物优化依赖资深化学家的经验,拖慢了药物开发进程。AI生成化学空间的快速扩张需要高效的可解释AI方法,这些工具能够发现超越人类直觉的洞见,补充专家知识并显著加速优化周期。为应对这些挑战,本文提出了一种新颖的双阶段原子归因值计算框架:包括基于结构池化的模型训练和基于子结构映射的原子归因值计算。该可解释框架量化任务特定的原子归因值,在芳香性/LogP/TPSA数据集上使用GCN时,原子归因准确性(计算值与真值的一致性)分别提升30%/20%/15%,Pearson相关系数达到0.93/0.63/0.88,超越了常用可解释方法仅能达到的0–0.3。此外,该方法对模型参数变化不敏感,对化合物结构变化提供相对稳定的预测结果。 核心结论 子结构池化+映射实现全局可解释的原子归因,显著提升与化学真值的一致性。 在芳香性、LogP、TPSA三任务上,TradePool的GCN原子归因Pearson相关0.93/0.63/0.88,F1、sparsity等指标全面优于GNNExplainer、KernelSHAP、Integrated Gradients、PGMExplainer,解释精度与稀疏性双优。 低频子结构筛除(出现次数<100)可抑制过拟合,保证权重的统计显著性。 对模型超参数和输入分子微扰不敏感,归因稳定性优于对照方法;但在GAT上效果一般,暴露了注意力权重与子结构加权的不匹配。 PubChem指纹提供任务无关的标准子结构集合,便于跨数据集、跨架构复用,部署与迁移成本低。 背景 图神经网络通过消息传递捕捉分子拓扑,在溶解度、毒性、反应性等性质预测上已成为主力。但多层聚合带来的“黑盒”问题削弱了可信度,尤其在药物优化环节,需要知道哪几个原子驱动了预测。 现有解释方法存在三大痛点:局部性强,难得到全局稳定的原子归因;与化学真值偏差大,Pearson相关常徘徊在0–0.3;计算代价高或对超参数敏感。子结构层面的解释更接近化学直觉,但GNN输入并未直接包含预定义子结构,如何把“可解释的子结构权重”映射回原子,成了瓶颈。 关键科学问题 如何在不牺牲预测精度的前提下,将GNN的决策过程转化为“子结构→原子”的可量化归因? 子结构集合应如何选择,既具普适性又能捕捉任务相关模式? 归因结果能否对模型参数、输入扰动保持稳定,从而在真实药物优化中可复用? 创新点 双阶段框架:训练时用PubChem指纹做子结构池化,解释时把子结构权重映射为原子归因。 全局归因:通过线性层权重直接量化子结构重要性,再按子结构-原子掩码汇总为原子级贡献。 稳健性设计:低频子结构剔除、权重聚合、多任务对比,提升对超参数和分子扰动的鲁棒性。 任务通用性:同一套指纹子结构跨芳香性、LogP、TPSA乃至药物临床分子数据集均可复用。 研究内容 方法详述 TradePool的核心思想是将子结构作为连接原子和分子性质的桥梁。传统GNN直接从原子嵌入池化到分子表示,丢失了化学家熟悉的官能团或子结构这一中间层信息。TradePool通过引入PubChem指纹定义的881个标准子结构,在训练时显式地学习每个子结构对预测的贡献权重,在解释时将这些权重映射回原子,从而实现全局一致、化学可解释的原子归因。 数据准备与清洗 为什么需要严格的数据清洗? 分子数据常存在SMILES表示不规范、含盐、带电荷等问题,这些会导致同一分子有多种表示形式,影响模型训练和归因评估的准确性。 研究使用RDKit 2022.09.5和MolVS 0.1.1进行标准化处理,包括SMILES规范化统一分子表示确保同一分子只有唯一的SMILES字符串,去盐处理移除分子中的无机盐(如$\ce{NaCl}$、$\ce{HCl}$)只保留有机部分,中和处理将带电荷的分子转为中性形式避免电荷状态影响特征计算,以及去重按分子骨架去除重复化合物防止数据泄漏。 中和应该存疑,应该是所选pH下的状态 清洗后的数据按8:1:1比例划分为训练集、验证集和测试集,这种划分确保模型在训练时不会接触测试集分子,从而真实评估泛化能力。 特征工程:从分子到图 原子特征(71维):每个原子用71维向量描述,包含11类信息 原子类型(43维,C、N、O、S等元素的one-hot编码)、度数(11维,原子连接的其他原子数量0-10+)、隐式价(7维,未显式表示的氢原子数) 电荷(原子的形式电荷如-1、0、+1)、芳香性(是否为芳香原子)、自由基电子(未配对电子数)、杂化类型(sp、sp²、sp³等)、连接氢数(显式连接的氢原子数)、手性中心(是否为手性中心)、手性类型(R/S构型) 键特征(12维)包含4类信息:键型(4维,单键、双键、三键、芳香键)、共轭性(是否参与共轭体系)、是否在环中(环状结构标识)、立体化学(E/Z构型或顺反异构) 这些特征由RDKit自动计算,详见Supporting Information Table S1。 子结构筛选:从881位到400+位 为什么要筛选子结构? PubChem指纹包含881个预定义子结构,但并非所有子结构都在数据集中频繁出现。低频子结构(出现次数<100)在统计上不显著,可能导致模型过拟合——模型会记住这些稀有模式而非学习真正的化学规律。 筛选策略:统计每个子结构在数据集中的出现频次,设定阈值为出现次数≥100次才保留,移除低频子结构以减少噪声和过拟合风险。 筛选结果(图1):芳香性任务保留416个子结构(移除465个),LogP任务保留513个子结构(移除368个),TPSA任务保留442个子结构(移除439个) 图1:三类任务的子结构出现频率热图 横轴:PubChem指纹子结构位;纵向颜色深浅:出现频率占比,深色表示更常见 截断低频(<100次)后,仍可覆盖大多数分子,避免稀疏噪声 筛选后仍能覆盖>90%的分子,说明被移除的子结构确实是稀有模式。图1的热图显示,保留的子结构在数据集中分布相对均匀,颜色深浅代表出现频率——深色表示高频子结构,浅色表示中频子结构。 数据集与标签构建 研究选择了三个具有明确原子归因真值的任务,这是评估XAI方法准确性的关键——只有存在可对照的真值,才能判断模型的解释是否可信。 芳香性数据集(Aromaticity) 为什么选择芳香性? 这是唯一具有客观真值的数据集,被多篇XAI论文用作基准测试。芳香性是分子的固有结构属性,不依赖于计算方法,一个原子是否芳香可以通过Hückel规则明确判定,因此原子归因的真值是确定的。 数据集来源:本文沿用Xiong等人构建的芳香性数据集,用于检验模型在原子层面的化学可解释性,标签为每个分子中芳香原子的数量。 LogP数据集(脂溶性) LogP的化学意义:LogP衡量分子的疏水性,是药物设计中的关键参数。疏水性高的分子更容易穿透细胞膜,但过高会导致溶解度差。 数据集来源:本文使用Wang等人整理的脂溶性数据集,主要来自PHYSPROP数据库与Hansch汇编数据集。 原子归因真值:以Crippen方法给出的原子级LogP贡献作为真值,用于评价连续归因值与真值的一致性。 额外外部集合:411个FDA批准药物与10个SAMPL6挑战分子被用作外部评估,用于检验不同方法的原子归因效果在真实药物结构上的表现。 Crippen原子贡献法是什么 TradePool把Crippen方法当作LogP任务的“原子归因真值”,这一步非常关键,因为它让“解释对不对”变成了可量化的问题。 方法来源:Wildman与Crippen在1999年提出一种原子类型分类体系,用原子贡献加和来预测分子的logP与摩尔折射率(MR)。 核心思想:先根据每个原子的局部化学环境把它分到某个原子类型,再把对应类型的贡献值相加得到全分子的logP。 计算形式:分子的logP可写作 \(\log P = \sum_{i=1}^{N} a_{t(i)}\) 其中,$t(i)$表示原子$i$所属的原子类型,$a_{t(i)}$是该类型的经验贡献系数,$N$是原子数。 为什么适合作为“真值”:它天然给出每个原子的数值贡献,可直接与XAI输出的连续归因值做Pearson相关比较。 RDKit里的实现:RDKit在rdkit.Chem.Crippen模块中提供MolLogP与MolMR,明确采用Wildman–Crippen的原子贡献方案;计算时还提供addHs选项,允许在需要时临时补氢参与贡献计算。实际结果会受到芳香性判定与是否显式加氢的影响,因此同一SMILES在不同标准化流程下可能出现轻微差异。 需要牢记的局限:Crippen是经验模型,主要面向中性小分子;它描述的是分子在辛醇与水相之间的分配倾向,不直接等同于带电体系的logD,也不显式建模溶剂化与构象效应。 参考:Wildman, S. A.; Crippen, G. M. Prediction of Physicochemical Parameters by Atomic Contributions. J. Chem. Inf. Comput. Sci. 1999, 39, 868–873. https://doi.org/10.1021/ci990307l TPSA数据集(拓扑极性表面积) 为什么TPSA重要? TPSA是药物类药性的关键指标,能够预测药物的溶解度、渗透性和药代动力学性质。一般认为,TPSA小于140 Ų的分子更容易口服吸收;极性表面积过大的分子难以穿透肠道上皮细胞,导致口服生物利用度降低。 数据与真值口径:TPSA本质上是一个基于分子拓扑的分子描述符,经典定义来源于Ertl等人的碎片贡献思想。本文将TPSA作为预测标签,并以碎片贡献法得到的原子级贡献作为归因真值,用于量化解释的正确性。 临床分子集(Drug-like Compounds) 为什么需要临床分子集? 前两个数据集虽然有真值但分子多样性有限,临床分子集包含真实的II期及以上候选药物,骨架复杂度更高,更能测试TradePool在实际药物优化场景中的可迁移性。 数据集来源:作者从ChEMBL数据库收集5800个分子量0到600、处于II期及以上临床阶段的小分子;按骨架结构划分训练、验证与测试集。该任务的训练标签与原子归因真值均由RDKit计算。 表1 数据集关键信息对比 | 数据集 | 总样本量 | 训练/验证/测试 | 任务类型 | 原子归因真值 | 数据来源 | 额外测试集 | | — | — | — | — | — | — | — | | 芳香性 | 3947 | 3157/395/395 | 回归(芳香原子数) | 芳香原子标签(芳香原子为1) | Xiong等构建 | - | | LogP | 16296 | 13036/1630/1630 | 回归(辛醇-水分配系数) | Crippen原子贡献 | PHYSPROP与Hansch汇编(Wang等整理) | 411个FDA药物,10个SAMPL6 | | TPSA | 5800 | 4700/550/550 | 回归(拓扑极性表面积) | 碎片贡献法原子贡献 | 文中未详述分子来源 | - | | 临床分子 | 5800 | 按骨架划分 | 由RDKit计算的分子性质 | RDKit计算 | ChEMBL(II期及以上) | - | 模型架构与训练策略 GNN编码器选择 研究实现了三种主流GNN架构,以验证TradePool的通用性: GCN(图卷积网络):每个原子聚合其邻居原子的特征,通过加权求和更新自身表示。GCN简单高效,适合捕捉局部拓扑结构,实现基于PyTorch和DGL-LifeSci 0.3.2。 GraphSAGE(图采样聚合):采样固定数量的邻居,使用LSTM或mean聚合器处理邻居特征。GraphSAGE可扩展到大规模图,聚合方式更灵活,超参数偏好LSTM或mean聚合器。 GAT(图注意力网络):为每个邻居分配注意力权重,动态调整不同邻居的重要性。GAT能够自适应地关注重要邻居,超参数采用4-8个注意力头,小或零dropout。 所有模型使用1-3层消息传递层,ReLU激活函数,隐藏维度在128-256之间。控制组使用传统的WeightedSumAndMax池化,TradePool组替换为子结构池化。 子结构池化机制 这是TradePool的核心创新。传统池化直接将所有原子嵌入求和或取最大值,丢失了子结构信息。TradePool的池化过程如下: 步骤1:构建子结构-原子掩码矩阵S。矩阵维度为$N \times P$,其中$N$是分子中的原子数,$P$是保留的子结构数(416/513/442),矩阵元素$S_{ij} = 1$表示原子$i$属于子结构$j$否则为0。计算方式使用RDKit的PubChem指纹生成函数,自动识别每个原子属于哪些子结构。 步骤2:创建子结构虚拟节点。对于每个子结构$j$创建一个虚拟节点$V_j$,虚拟节点的特征等于所有属于该子结构的原子嵌入之和:\(V_j = \sum_{i: S_{ij}=1} h_i\)其中$h_i$是原子$i$经过消息传递后的嵌入向量。如果分子不包含子结构$j$,则$V_j$为全零向量。 步骤3:展平与预测。将所有子结构虚拟节点展平为一维向量 $[V_1, V_2, …, V_P]$,输入到单层线性层$\hat{y} = W \cdot [V_1, V_2, …, V_P] + b$,输出预测标签(回归任务输出标量,分类任务输出类别概率) 为什么这样设计有效? 子结构池化强制模型通过子结构这一中间层进行预测,使得线性层的权重$W$直接对应每个子结构的重要性。这种设计天然地将可解释性嵌入模型架构,而非事后添加。 图2:TradePool双阶段工作流 (A) 总览:左侧训练阶段输入分子图与子结构掩码,右侧解释阶段输出原子归因 (B) 子结构池化:同一子结构内的原子特征求和形成虚拟节点;未包含该子结构则为零向量 (C) 归因映射:线性层得到子结构权重,按掩码回分到原子,权重累加得到原子归因热图 训练超参数与优化策略 优化器与学习率:使用Adam优化器自适应调整每个参数的学习率,学习率通过贝叶斯优化在验证集上搜索最优值,典型范围为1e-4到1e-3。批大小设定为128平衡内存占用与梯度稳定性,最大训练轮次为200,早停策略监控验证集MAE(回归)或准确率(分类),连续10轮无改善则停止。权重初始化使用Xavier初始化,确保每层输出方差一致,避免梯度消失或爆炸。 训练稳定性技巧:采用冻结策略,训练后10%的轮次仅微调线性层保持图编码部分冻结,目的是降低梯度震荡确保子结构权重稳定可解释。必要时引入L2正则化抑制极端权重,防止单一高频子结构独占权重。 超参数搜索结果(SI Table S2):隐藏维度为128-256,层数为2-3层,GraphSAGE偏好LSTM或mean聚合,GAT采用4-8个注意力头配合小或零dropout。TradePool与控制组使用相同深度,主要区别在池化方式。 原子归因计算 训练完成后,如何从子结构权重得到原子归因?这是TradePool的第二阶段——解释阶段。 提取子结构权重 不同GNN架构的权重提取方式不同,因为它们的聚合机制不同: GCN:线性层权重矩阵$W$的每一列对应一个子结构,子结构$j$的归因值等于该列所有元素之和$A_j = \sum_k W_{kj}$,原理是GCN的聚合是简单求和,权重的和反映了子结构的总贡献。 GraphSAGE和GAT:子结构$j$的归因值等于该列所有元素的L1范数$A_j = \sum_k W_{kj} $,原理是这些模型的聚合更复杂(LSTM或注意力),权重可能有正负,取绝对值后求和更稳定。 映射到原子归因 有了每个子结构的归因值$A_j$,如何得到每个原子的归因值? 映射公式为: \(\text{原子}i\text{的归因} = \sum_{j: S_{ij}=1} A_j\) 通俗解释:查找所有包含原子$i$的子结构(即$S_{ij}=1$的子结构),将这些子结构的归因值累加,累加结果即为原子$i$的归因值。化学直觉:同一子结构内的原子获得相同的基础贡献(因为它们都属于该子结构),处于多个子结构交叉位置的原子累积多重贡献(例如苯环上的碳既属于“芳香环”子结构,也属于“C=C”子结构),这种累加方式与化学家的思维一致——一个原子的重要性取决于它参与了哪些官能团。 呃,其实也可以在搞原子对子结构贡献的权重的,就有点复杂了 “正归因原子”如何定义:阈值与二值化 很多指标(Accuracy、Recall、F1、Sparsity)要求先把连续归因值变成二分类标签。论文对不同方法的二值化规则如下: 传统XAI方法:若原子归因值大于0,则标记为正归因;否则为负归因。 TradePool方法:由于原子归因值来自子结构权重累加,作者不直接使用0作为阈值,而是在训练集上计算一个任务级阈值:对训练集中每个分子,记录其原子归因最大值与最小值;对所有分子的最大值与最小值分别取均值;再取这两个均值的平均作为阈值。验证集与测试集沿用训练集得到的阈值。 通俗解释:这个阈值更像是在训练集的归因值动态范围里取一个平均中线,用它来区分相对更重要与相对不重要的原子。它并不强制每个分子都选出固定比例的原子,只是实际结果常落在中等稀疏度区间。 同时,论文也把真值归因二值化用于分类指标计算: 芳香性任务:芳香原子真值标签为1,其他为0。 LogP与TPSA任务:真值原子贡献大于0标为1,否则为0。 评估指标体系 研究采用Wang等人提出的八项XAI评估指标,全面衡量归因质量: 准确性指标: Accuracy:分类任务,正确识别正/负归因原子的比例 F1-score:精确率和召回率的调和平均,平衡误报和漏报 Pearson相关:预测归因值与真值的线性相关性(-1到1,越接近1越好) 稀疏性指标: Sparsity:被标记为正归因的原子比例。理想的解释应该聚焦于少数关键原子,而非高亮整个分子 Recall:真正的正归因原子中被正确识别的比例 稳定性指标: Fidelity:移除正归因原子后,预测值下降的幅度。下降越多,说明这些原子确实重要 Infidelity:移除负归因原子后,预测值上升的幅度。上升越多,说明这些原子确实有负贡献 Stability:对模型参数微调(如改变随机种子)后,归因结果的一致性 Sensitivity:对输入分子微扰(如添加甲基)后,归因结果的稳定性 呃,不一定非得有下降和上升的幅度很突出的原子吧,比如烷烃,都一样? 这些指标从不同角度评估XAI方法:准确性衡量解释是否正确,稀疏性衡量解释是否简洁,稳定性衡量解释是否可靠。只有在所有维度都表现优秀,才能称为真正好的XAI方法。 总结:TradePool双阶段框架 graph TB subgraph S1["阶段1 训练:子结构池化"] A1["分子图输入<br/>原子71维+键12维"] --> B1["消息传递<br/>GCN / GraphSAGE / GAT"] B1 --> C1["子结构池化<br/>PubChem掩码求和虚拟节点"] C1 end subgraph S2["阶段2 解释:权重映射"] D1["线性层预测标签<br/>同时产生子结构权重"] --> E1["选取子结构权重<br/>不同模型取和或L1范数"] E1 --> F1["按掩码回分原子<br/>归因=所有含该原子的子结构权重之和"] F1 --> G1["输出原子归因热图<br/>稳定、可化学解释"] end S1 --> S2 图2详细展示了TradePool的核心工作流程。TradePool的创新在于将子结构作为可解释的中间层,连接原子级输入和分子级预测。 训练阶段(图2A左侧):分子图经过消息传递层后,不是直接进行全局池化,而是根据预先计算的PubChem指纹子结构掩码,为每个子结构创建虚拟节点。这一步骤(图2B)通过将属于同一子结构的所有原子嵌入求和来实现——如果某个分子不包含某个子结构,则对应的虚拟节点为零向量。这些子结构表示随后被展平并输入到线性层以预测最终标签。 解释阶段(图2A右侧):TradePool的优势得以充分体现。由于线性层的权重直接对应于每个子结构对预测的贡献,研究者可以提取这些权重作为子结构归因值。图2C展示了如何将子结构归因映射回原子:对于每个原子,查找所有包含该原子的子结构,将这些子结构的权重累加,即得到该原子的归因值。这种从子结构到原子的映射策略确保了归因的全局一致性——同一子结构内的原子获得相同的基础贡献,而处于多个子结构交叉位置的原子则累积多重贡献,这与化学直觉高度一致。 结果与分析 主任务预测性能:并未牺牲预测精度 在GCN上,TradePool的原子归因Pearson相关:芳香性0.93,LogP 0.63,TPSA 0.88;常见解释方法多在0–0.30之间。 F1与sparsity均优于GNNExplainer、KernelSHAP、Integrated Gradients、PGMExplainer,说明归因更聚焦、冗余更少,解释“准而简”。 GAT上表现一般,源于注意力权重为标量,难与子结构权重对齐,提示池化假设需与注意力机制协同。 预测精度与对照组相当或更优:TPSA任务GCN的MAE 1.157、RMSE 1.569(对照3.367、4.846),LogP任务GCN的MAE 0.299与对照0.296持平,显示可解释性增强未牺牲主任务性能(SI Tables S3–S5)。 主任务预测性能是可解释性的基础。图3展示了TradePool在三个任务上的预测精度散点图,每行对应一个任务(芳香性、LogP、TPSA),每列对应一个GNN架构(GCN、GraphSAGE、GAT)。 从图中可以看到,所有模型在三个任务上都取得了较高的$R^2$值,数据点紧密分布在对角线附近,表明预测值与真实值高度吻合。值得注意的是,TradePool引入的子结构池化机制并未牺牲预测精度——在大多数情况下,TradePool的$R^2$与使用传统WeightedSumAndMax池化的对照组相当,甚至在TPSA任务上表现更优。这证明了子结构池化不仅提升了可解释性,还通过结构化的中间表示增强了模型对任务相关模式的学习能力。 图3:三任务的真值与预测散点 行:芳香性、LogP、TPSA;列:GCN、GraphSAGE、GAT 颜色区分训练/验证/测试;对角越集中表明拟合越好,右上角图例标示$R^2$ 原子归因精度全面领先 在保证预测性能的前提下,TradePool在原子归因质量上实现了显著突破。图4对比了TradePool与四种主流XAI方法(GNNExplainer、KernelSHAP、Integrated Gradients、PGMExplainer)在原子归因准确率和Pearson相关性上的表现。 图4A显示,TradePool在芳香性、LogP和TPSA三个任务上,特别是在GCN架构上,原子归因准确率均达到最高,芳香性任务甚至接近0.9的准确率。图4B的Pearson相关更是揭示了TradePool的优势:在GCN上,TradePool的相关系数达到0.93(芳香性)、0.63(LogP)和0.88(TPSA),而其他方法大多徘徊在0–0.3之间,部分方法甚至出现负相关。这种量级的提升表明,TradePool生成的原子归因不仅在统计上与化学真值一致,而且能够准确捕捉到原子对分子性质的实际贡献。 图4:原子归因总体表现对比 (A) 原子归因准确率:TradePool在三任务、特别是GCN上最高 (B) Pearson相关:TradePool显著领先,其他方法多在0–0.3之间 图5想回答的问题是:如果一个方法把很多原子都判成关键原子,它当然容易拿到高召回,但这不一定是好解释。作者强调要同时看F1、召回与稀疏性,避免靠把整分子都高亮来“刷指标”。 结论1:TradePool的召回不是靠过度归因换来的。论文指出,TradePool在芳香性与LogP任务的F1与召回都表现良好,更重要的是稀疏性维持在0.4–0.5,意味着大约只有40%–50%的原子被标为正归因,解释更聚焦、信息密度更高。 结论2:KernelSHAP与PGMExplainer存在明显的过度归因倾向。它们在部分任务上召回较高,但对应稀疏性很低,说明方法倾向把接近90%的原子都判为正归因,从而抬高召回。作者认为稀疏性过高或过低都意味着解释存在缺陷:太低会导致解释冗余,太高又容易漏掉关键结构片段。 图5:F1、召回与稀疏性对比 (A) F1-score:TradePool在LogP/TPSA上优势明显。 (B) 召回率:KernelSHAP与PGMExplainer高召回但伴随过多正归因。 (C) 稀疏性:TradePool保持0.4–0.5的稀疏度,解释更集中。 为了直观展示不同XAI方法的归因质量,研究团队随机选取了三个任务(芳香性、LogP、TPSA)测试集中的分子,使用GCN模型生成原子归因热图并进行对比。图11中,绿色高亮表示被预测为正归因的原子,每一行对应一个不同的任务。 从可视化结果可以清晰看到,TradePool的原子归因与化学真值高度吻合:在芳香性任务中,TradePool准确高亮了芳香环上的原子;在LogP任务中,疏水性基团(如芳环、烷基链)被正确识别;在TPSA任务中,含氧、含氮的极性原子得到强调。相比之下,GNNExplainer、KernelSHAP、Integrated Gradients和PGMExplainer等方法存在明显的误高亮问题——它们倾向于高亮更多的原子,包括一些与任务无关的位点,导致归因稀疏性降低、解释冗余增加。 这种可视化案例验证了前面定量指标的结论:TradePool不仅在Pearson相关、F1等数值指标上优于对照方法,在实际化学结构解释的视觉一致性上也表现更佳,更符合化学家的直觉判断。 图11:不同XAI方法在GCN模型上的原子归因可视化对比 每行对应芳香性、LogP和TPSA三个任务之一。 绿色高亮表示被预测为正归因的原子。 TradePool的高亮区域与真值最吻合,误高亮最少;其他方法存在明显的过度归因或归因偏差。 子结构化学合理性 提取各任务权重前十的SMARTS子结构(表3),与化学常识一致: 芳香性任务40%含芳香键; LogP任务突出甲基、芳环、卤素; TPSA任务高频出现含氧、含氮片段O−H、N−O、S(=O)(=O)。 子结构权重跨随机种子保持一致(SI Table S6),N=N、C=S、N−S、O(:C)(:C)等始终位列前十,支持归因的可重复性。 在SMARTS里,: 表示芳香键(aromatic bond),:C 表示芳香碳原子(aromatic carbon) ~ 表示任意键(any bond):不限定是单键、双键、三键还是芳香键,只要两原子之间“有键”就匹配。 表3 三个任务权重最高的10个SMARTS子结构 芳香性 LogP TPSA N#N C−I ≥1 O C=S ≥1 Br O(:C)(:C) N−S N#N N−H C(∼N)(:N) ≥1 Cl ≥1 P N−C:O:C ≥1 S O−H C(:N)(:N) C:C−N−C:C N#N ≥1 Cl C−Br ≥1 N O(:C)(:C) C(∼H)(∼H)(∼H) C#N N−N C(∼F)(∼F) S(=O)(=O) C#N ≥1 F N=O 稳定性与鲁棒性 对模型参数微调或输入分子小幅扰动,TradePool的fidelity/infidelity方差最低,归因热图变化最小。 归因稳定性的原因:子结构集合固定、权重全局学习、低频子结构滤除减少噪声。 对指纹掩码随机置零10%或对分子增加单键旋转等扰动,原子归因排名的Spearman相关仍高于0.85,而对照方法掉到0.5以下,说明结构微扰下解释更稳。 在超参数网格搜索(学习率、隐藏维度、层数)中,TradePool的Pearson相关标准差小于0.03,显著优于对照,超参敏感性低。 Fidelity/Infidelity统计(SI Tables S7–S8):LogP任务TradePool的fidelity_mean=3.38、infidelity_mean=0.69,明显优于其他方法;TPSA任务TradePool保持正向fidelity 36.18,而KernelSHAP虽fidelity高但infidelity为负且方差大,说明TradePool稳定性更高。 图6对应论文的稳定性实验:作者在芳香性测试集里对369个含苯环的分子做结构微扰,在苯环上添加1–2个甲基;不含苯环的分子不做修改。随后用同一个预训练模型分别对加甲基前后分子计算原子归因,并比较每个分子的原子归因准确率变化量。 结论:TradePool对结构小改动更稳。论文报告所有方法的变化总体接近0,但TradePool的变化显著更小,说明当分子发生轻微修饰时,TradePool给出的关键原子集合更不容易漂移。对于药物优化而言,这意味着解释可以跨相邻类似物复用,降低“每做一次修饰就要重新理解解释”的成本。 为什么这能叫稳定:芳香性任务的真值关键原子主要是芳香环原子,给苯环加甲基不会改变原来的芳香原子标签。理想的解释应继续高亮芳香环,而不是被新增甲基带跑偏。 图6:小扰动下的稳定性 (A) 在分子上随机添加1–2个甲基的示例。 (B) 各方法扰动前后原子归因值变化,TradePool波动最小,说明对结构微扰不敏感。 图7是图6的可视化证据:同一批分子在加甲基前后的归因热图对比。这里的绿色代表被方法预测为正归因的原子,也就是它认为的关键原子。在芳香性任务里,这些原子理想情况下应与真实的芳香原子位置一致。 结论1:TradePool的高亮区域更贴近化学真值且更一致。加甲基前后,TradePool主要持续高亮芳香环原子,新增甲基不会导致模型把大量非芳香原子误判为关键。 结论2:部分对照方法会把几乎所有原子都判为芳香。论文特别指出KernelSHAP与PGMExplainer会把所有原子都分类为芳香原子,这会造成“看似召回很高、实际毫无区分度”的解释,和图5中稀疏性异常的问题相呼应。 把图5–7连起来读:TradePool不仅在指标上避免过度归因,还能在结构微扰后保持解释形状;而某些方法的高召回来自过度归因,导致热图失去化学可读性。 图7:扰动前后的原子热图对比 绿色高亮:被判定为正归因的原子;每两行对应同一分子扰动前后。 TradePool在扰动后保持高亮区域一致性,对照方法出现更多误高亮。 图8:跨随机种子的敏感性。五个随机种子训练的GCN模型在芳香性测试集的原子归因准确率箱线图;TradePool方差最小,鲁棒性最佳。 图9:不同随机种子下的原子热图 每行对应一个随机种子训练的模型,绿色为正归因原子。 TradePool跨种子保持高亮模式一致,对照方法高亮位置漂移更大。 Fidelity的实验验证是通过移除原子来测试归因质量:将模型预测为正归因的原子张量置零后重新预测,若预测值显著下降,说明这些原子确实对预测有正贡献;反之,移除负归因原子应使预测上升。 图10展示了在芳香性任务上,TradePool移除正归因原子后预测值下降最多,移除负归因原子后预测值上升最稳定,证明其归因方向与化学机制高度一致。相比之下,其他方法如GNNExplainer和Integrated Gradients移除原子后预测值变化较小或方向不一致,表明归因质量不佳。 纵轴:模型对芳香性任务的输出(预测的芳香原子数量),数值越大表示预测的芳香原子越多。 横轴(1–4)与图例一一对应:1为真值,2为完整分子预测,3为仅保留关键原子后的预测,4为移除关键原子后的预测。 若一个方法真的找到了关键原子,那么黄色箱线图应明显低于绿色(拿走关键后预测下降),而红色应仍接近绿色(只看关键也能维持预测)。 图例解释:蓝色为真值分布;绿色为完整分子预测;红色表示把非关键原子特征置零后的预测;黄色表示把关键原子特征置零后的预测。 图10:移除正/负归因原子后的预测值变化 在芳香性任务上,移除不同XAI方法预测的正/负归因原子后,模型的预测值变化。 TradePool移除正归因原子后预测值下降幅度最大(KernelSHAP也还行?),移除负归因原子后预测值上升幅度最小且稳定,验证其归因方向正确。其他方法移除原子后预测值变化较小或方向不一致,说明归因与模型决策机制不匹配。 版面结构:每一列对应一种解释方法(Random、TradePool、GNNExplainer、KernelSHAP、IG、PGMExplainer),每一行对应一种GNN架构(GCN、GraphSAGE、GAT)。箱线图显示分布范围,箱体中线为中位数,三角形为均值。 临床分子集迁移 在5800个II期及以上候选药物上,TradePool在GCN与GAT上均优于传统基线,显示其对真实药物骨架的可迁移性。 典型案例:含卤代芳环的口服候选物,TradePool高亮芳环与卤素原子,与疏水性主导的LogP真值一致;对照方法偏高亮杂原子,解释偏差较大,显示对实际药物骨架的解释可靠性。 在LogP与TPSA任务中,去除正归因原子会导致TradePool预测下降幅度最大,去除负归因原子则上升最小(SI Figures S1–S2),进一步验证其归因方向符合化学机制。 需要注意的是:正文与Supporting Information未给出该临床分子集的完整数值表,仅给出数据集构建方式与文字性结论描述。 结果逻辑图(方法—结果—局限) graph TB subgraph S0["问题与设计"] Q1("难以获得稳定原子归因") --> Q2("采用指纹子结构做全局池化") end subgraph S1["实验管线"] D0("数据清洗与指纹筛选") --> D1("GCN/GraphSAGE/GAT训练") D1 --> D2("线性层子结构权重") D2 --> D3("掩码回分原子归因") end subgraph S2["核心结果"] R1("Pearson相关0.93/0.63/0.88") --> R2("F1与sparsity领先基线") R2 --> R3("扰动下归因稳定性最佳") R1 --> R4("权重前十子结构符合化学直觉") end subgraph S3["局限与改进"] L1("与GAT注意力不匹配") --> L2("计划加入可学习子结构生成") L3("小数据高复杂任务R^2偏低") --> L2 end Q2 --> D0 D3 --> R1 D3 --> R3 R3 --> L1 方法局限与改进方向 对注意力模型支持不足:子结构等权假设与GAT的原生注意力冲突。 数据集较小(<3000)或任务复杂时,$R^2$偏低,子结构权重难以学到任务相关性。 未来计划:在训练中加入“生成-筛选”子结构模块(类似GAN),替换低权重子结构,提升任务相关性与多样性。 化学与工程解读 化学角度:子结构权重凸显芳香键、卤素、含氧氮片段,与芳香性、疏水性、极性表面积的主导因素一致,提升了模型的化学可信度。 工程角度:使用标准指纹可避免任务特定规则,部署时只需计算指纹与权重矩阵,无需逐分子重新训练,适合大规模虚拟筛选。 Q&A Q1:为什么用PubChem指纹而不是ECFP或规则切分? A1:PubChem指纹是公开字典,881位覆盖常见官能团,跨分子可比;数量适中,便于全局权重学习;规则切分在多数分子下碎片数<10,统计显著性不足。 Q2:子结构权重如何转成原子归因? A2:训练后从线性层取每个子结构的权重(GCN取和,GraphSAGE/GAT取L1范数),再用子结构-原子掩码,将包含该原子的所有子结构权重相加,即为该原子的归因值。 Q3:为什么对参数和分子扰动更稳? A3:归因依赖全局训练得到的固定子结构权重,而非逐样本优化;子结构数量大、权重聚合降低单一掩码变化带来的波动;低频子结构被剔除减少噪声。 Q4:数据清洗如何保证标签一致性? A4:使用RDKit与MolVS标准化SMILES、去盐和中和,重复分子按骨架去重;标签计算遵循Crippen原子贡献或拓扑表面积分拆,保证训练与真值口径一致。 关键结论与批判性总结 潜在影响:为分子GNN提供全局、量化的原子归因路径,能直接指导先导优化与毒性定位,降低对专家经验的依赖。 局限性:与注意力类模型存在机制不匹配;小数据、高复杂任务下权重难学;对子结构词表的覆盖度仍依赖预定义指纹。 未来方向:引入可学习的子结构生成与淘汰机制;探索与GAT兼容的子结构加权方式;将方法拓展到蛋白-配体复合物、材料晶格等更大图结构。 小编锐评: 做可解释性分析的一种尝试了。我的体会是,如果更贴近人类语言,那还得是基团,但到底谁贡献多,会不会有相关,本身就是有点复杂的,case by case的解释是避免不了的。现在这样有解释已经不错了。 做可解释性分析可以水这么多图,学到了
Machine Learning & AI
· 2026-01-11
黑箱的透明化:通过原子敏感性分析实现可解释的pKa预测模型
“黑箱”的透明化:BCL-XpKa通过原子敏感性分析实现可解释的pKa预测模型 Title: Interpretable Deep-Learning pKa Prediction for Small Molecule Drugs via Atomic Sensitivity Analysis Authors: Joseph DeCorte,* Benjamin Brown, Rathmell Jeffrey, and Jens Meiler https://doi.org/10.1021/acs.jcim.4c01472 Cite This: J. Chem. Inf. Model. 2025, 65, 101-113 一、 论文整体概览 1. 摘要、背景与科学问题 摘要翻译 机器学习(ML)模型如今在预测药物研发所必需的性质方面扮演着至关重要的角色,例如药物的对数尺度酸解离常数(pKa)。尽管近期在架构上取得了进展,但由于缺乏基准真实数据,这些模型在面对新化合物时常常泛化能力不佳。此外,这些模型也缺乏可解释性。为此,通过精心设计的分子嵌入,可以通过观察模型对输入分子进行原子扰动后的响应,来获取化学结构中的原子级分辨率信息。在此,我们提出了BCL-XpKa,一个基于深度神经网络(DNN)的多任务分类器,用于pKa预测,它通过Mol2D描述符来编码局部原子环境。BCL-XpKa为每个分子输出一个离散分布,该分布存储了pKa预测值以及模型对该分子的不确定性。BCL-XpKa能很好地泛化到新的小分子上,其性能与现代ML pKa预测器相当,在泛化任务中优于多个模型,并能准确模拟常见分子修饰对分子可电离性的影响。然后,我们通过原子敏感性分析(ASA)利用BCL-XpKa的精细描述符集和以分布为中心的输出,该分析无需重新训练模型即可将分子的预测pKa值分解为其各自的原子贡献。ASA揭示了BCL-XpKa已经隐式地学习到了关于分子亚结构的高分辨率信息。我们进一步通过在93.2%的复杂小分子酸和87.8%的碱中识别电离位点,展示了ASA在为蛋白质-配体对接准备结构方面的效用。最后,我们应用带有BCL-XpKa的ASA方法,识别并优化了一款最近发表的KRAS降解PROTAC的物理化学缺陷。 背景 在计算辅助药物研发领域,准确预测化合物在体内的行为(如生物利用度、溶解度等)对于节约研发时间和成本至关重要。其中,分子的酸解离常数(pKa)是一个决定其在生理pH下电离状态的关键物理化学性质,深刻影响着药物的吸收、分布、代谢、排泄和毒性(ADMET)。 传统的预测方法中,量子力学(QM)计算能够提供与实验相当的精度,但其巨大的计算成本使其无法应用于药物发现早期阶段对数以亿计化合物的虚拟高通量筛选(vHTS)。因此,机器学习(ML)方法,特别是定量结构-活性/性质关系(QSAR/QSPR)模型,因其极高的预测速度而成为主流。这些模型通过分子指纹或图神经网络(GNNs)等方式将化学结构转化为数学表示,并学习结构与性质之间的关系。 本文解决的科学问题总结 尽管ML方法取得了巨大成功,但仍面临两大核心挑战,这也是本文着力解决的科学问题: 性能与泛化问题:现有的ML模型大多依赖于数量有限的高质量实验数据进行训练,这常常导致模型在面对训练集中未见过的、新颖的化学骨架时泛化能力差,容易过拟合。 可解释性问题:大多数先进的ML模型(尤其是深度学习模型)如同一个“黑箱”,我们很难理解模型是基于分子的哪些具体结构特征做出某一特定预测的。这种可解释性的缺乏阻碍了我们对模型预测结果的信任,也使得我们难以从模型的“智慧”中获得化学洞见来指导后续的药物设计。 本文旨在通过创新的模型架构(BCL-XpKa)和新颖的可解释性分析方法(ASA)来同时应对这两个挑战。 mindmap root(可解释性pKa预测分析思路) )为可解释性服务的模型架构( ::icon(fa fa-cogs) **多任务分类(MTC)架构** **核心思路**<br/>预测pKa落入离散区间的概率分布 **关键产出**<br/>输出**概率分布**而非单个值<br/>为ASA提供可比较的分布 **附加价值**<br/>分布的标准差可作为**模型不确定度** **局部原子描述符(Mol2D)** **核心思路**<br/>仅编码原子及其一阶邻居<br/>使模型对单原子扰动更敏感 **关键特性**<br/>**可逆性**:描述符可直接映射回化学亚结构<br/>是实现归因分析的基础 )原子敏感性分析(ASA)( ::icon(fa fa-atom) **核心方法:基于扰动的归因** **扰动方式**<br/>将目标杂原子替换为同构的碳原子 **差异量化**<br/>用**KL散度**衡量扰动前后<br/>模型输出的pKa概率分布差异 **分数计算**<br/>通过指数函数放大KL散度<br/>得到最终ASA分数 **应用一:识别关键电离位点** **做法**<br/>寻找分子中ASA分数最高的原子 **结论**<br/>高分原子大概率是主要质子化或去质子化位点 **价值**<br/>快速标注分子质子态<br/>用于对接或MD模拟前的结构准备 **应用二:洞察模型学习到的化学知识** **做法**<br/>分析特定亚结构在不同化学环境下的ASA分数变化 **化学原理验证**<br/>**诱导效应**:邻近吸电子或给电子基团<br/>会相应降低或提高碱性氮的ASA分数 **上下文理解**<br/>模型能区分局部环境相同但整体有别的基团<br/>(如咪唑 vs. 吲哚) **研究启发**<br/>验证模型是否学到真实的化学规则<br/>增加对“黑箱”预测的信任 **应用三:指导先导化合物优化** **完整工作流**<br/>1. **定位缺陷**:用ASA找到导致不良pKa的原子<br/>2. **提出方案**:对高分原子进行生物电子等排替换<br/>3. **快速验证**:用模型预测新分子的pKa<br/>4. **结构确认**:用对接等方法确认活性 2. BCL-XpKa 模型简介 作者首先构建了一个名为 BCL-XpKa 的pKa预测模型,其核心是一个多层感知机(MLP)。该模型的设计巧妙,集成了几个关键特性: 分类而非回归(Multitask Classification, MTC):不同于传统模型直接预测一个连续的pKa值,BCL-XpKa将pKa范围划分为多个离散的“桶”(bins),并预测分子的pKa值落入每个“桶”的概率。最终的pKa值是这个概率分布的期望值。这种做法的好处是: 可以直接从输出分布的标准差中读出模型对预测的不确定度。 通过识别模型在哪些分子上表现出“高不确定性”或“高置信度但高错误率”,可以指导训练数据的优化。 在性能上与回归模型相当,甚至略优。 双模型架构:为了处理既有酸性基团又有碱性基团的复杂分子,作者分别训练了 BCL-XpKaAcid 和 BCL-XpKaBase 两个模型,用于分别预测一个分子中酸性最强和碱性最强的pKa值。 图1:BCL-XpKa的架构评估 (A) BCL-XpKa 使用独立的模型来预测分子的酸性和碱性pKa值。它使用 Mol2D 局部原子环境描述符来嵌入分子,然后使用一个多层感知机(MLP)来对pKa值所属的1-pKa单位区间进行分类。区间边缘交替包含和不包含端点。极值区间(pKa≤0, pKa>12)在其无界的一侧是开放的。 (B) 用于pKa预测的多任务分类误差随“桶”尺寸的增加而变化。小的“桶”允许更高的精度,但每个桶的数据更少;而大的“桶”精度较低,但每个桶的数据更多。 (C) BCL-XpKa与使用相同分子描述符和训练集训练的最佳性能回归架构在两个外部测试集上的性能对比。 (D) “留下一类” (leave-class-out, LCO) 方法,其中一种分子亚结构被从模型训练中移除,并在之后用作结构新颖的测试集。 (E) 模型误差由LCO亚结构和描述符类型决定。 (F) LCO亚结构的误差与包含该亚结构的TS-Acid或TS-Base分子数量的关系。 3. Mol2D 描述符为何对 ASA 至关重要? BCL-XpKa模型选择使用 Mol2D 描述符 而非更复杂的GNN,这是实现原子级别可解释性(ASA)的基石。参考其原始论文 BCL::Mol2D—a robust atom environment descriptor for QSAR modeling and lead optimization,Mol2D 的核心优势在于其设计上的简洁性与可逆性。 核心定义:Mol2D的核心是原子环境(Atom Environment, AE)。一个 AE 是以某个原子为中心,包含其周围一定化学键距离内的原子及其成键信息。BCL-XpKa使用的是 height=1 的AE,这意味着它只考虑中心原子和与它直接相连的邻居原子。 与传统指纹的关键区别: 计数而非存在与否:传统指纹(如 Molprint2D)通常是二进制的,只记录某种AE是否“存在”。而 BCL::Mol2D 是一个计数向量,它记录了分子中每种特定AE出现的次数。这提供了更丰富的信息,例如可以区分五元环和六元环。 细粒度的原子类型:Mol2D 不仅考虑元素类型,还考虑了原子的杂化状态/轨道构型(’Atom type’ 编码),这使得它能够区分同样是氮原子,但在不同化学环境下的细微差别。 通用 AE 库:BCL::Mol2D 的描述符向量的每一个维度都对应一个从大型化合物库(超过90万个类药分子)中预先构建好的“通用AE库”中的特定AE。这意味着描述符的索引是固定的,任何分子都可以被映射到这个统一的向量空间中。 可逆性(Reversibility)——实现ASA的关键: 这是 Mol2D 最重要的特性。由于描述符向量的每个索引都唯一地、固定地对应着一个具体的化学亚结构(即一个AE),我们可以从描述符向量反推回它所代表的化学结构。 这种清晰的“描述符-结构”对应关系,使得当我们扰动一个原子时,我们能精确知道是哪些维度的描述符发生了变化。这为衡量模型对特定原子变化的敏感度提供了直接、无歧义的途径。 相比之下,许多复杂的GNN模型其内部表示(节点嵌入)是经过多轮信息传递后高度抽象化的向量,难以直接映射回具体的、独立的原子或化学键贡献,从而使原子级别的归因分析变得非常困难。 4. BCL-XpKa 模型性能与表现总结 BCL-XpKa模型尽管采用了相对简单的多层感知机(MLP)架构,但在多个基准测试中展现了极具竞争力的性能。 与主流预测器的性能对比:在多个外部标准测试集(如Novartis、SAMPL6-8)上,BCL-XpKa的平均绝对误差(MAE)与包括ChemAxon、QupKake以及基于GNN的MolGpKa和Uni-pKa等在内的多种现代pKa预测器不相上下。例如,在Novartis-Acid测试集上,其MAE为0.79。 优秀的泛化能力:模型的核心优势在于其对新化学骨架的泛化能力。在“留下一类”(Leave-Class-Out, LCO)的交叉验证中,模型需要预测从未在训练集中见过的、特定化学亚结构分子的pKa。结果显示,使用Mol2D描述符的BCL-XpKa显著优于使用传统MACCS和Morgan指纹的同等模型,平均MAE分别为1.1(BCL-XpKa)、1.46(MACCS)和1.20(MFP2)。 准确捕捉化学趋势:模型不仅能预测绝对pKa值,更能准确地再现微小化学修饰所引起的pKa变化趋势。例如,在包含71对仅有细微结构差异的分子测试中,BCL-XpKa能够正确预测pKa变化方向的比例高达81.7%。这对于指导药物化学中的先导化合物优化尤为重要。 数据策略的有效性:该研究还表明,尽管模型主要在预测数据(来自ChEMBL)上进行训练,但其性能全面优于仅使用少量实验数据训练的同等模型(BCL-MLP-MTC-EO),证实了在大规模预测数据基础上进行训练策略的有效性。 二、 原子敏感性分析(ASA)方法细节与应用 这部分是该研究的核心。原子敏感性分析(Atomic Sensitivity Analysis, ASA) 的设计初衷是:在不重新训练模型的情况下,将模型对整个分子的pKa预测值“分解”到每个原子上,从而理解哪个原子或基团对最终的预测贡献最大。 1. ASA的核心原理 ASA的核心思想是“基于扰动的敏感性分析”。它通过系统性地、有物理意义地扰动分子中的每一个原子,并观察模型预测结果的变化剧烈程度,来判断该原子对原始预测的重要性。如果对某个原子的微小改动导致了模型预测结果的巨大变化,那么这个原子就被认为是“敏感的”或“重要的”。 2. ASA的具体实施步骤 graph TD subgraph "ASA 核心流程(针对单个原子)" A["**第1步:父本预测**<br/>将原始分子输入BCL-XpKa<br/>获得pKa概率分布 P_parent"] --> B; B["**第2步:原子扰动**<br/>将分子中的杂原子 a<br/>替换为同构的碳原子<br/>生成'扰动分子'"] --> C; C["**第3步:扰动预测**<br/>将'扰动分子'输入模型<br/>获得新的pKa概率分布 P_perturbed"] --> D; D["**第4步:量化差异**<br/>计算两个分布的差异<br/>使用Kullback-Leibler (KL)散度<br/>D_KL(P_perturbed || P_parent)"] --> E; E["**第5步:计算ASA分数**<br/>通过经验公式放大信号<br/>ASA(a) = exp[S * D_KL] - 1<br/>得到原子 a 的敏感性分数"] end style A fill:#e3f2fd,stroke:#1e88e5,stroke-width:2px style B fill:#fff3e0,stroke:#ef6c00,stroke-width:2px style C fill:#e3f2fd,stroke:#1e88e5,stroke-width:2px style D fill:#e8f5e9,stroke:#2e7d32,stroke-width:2px style E fill:#f3e5f5,stroke:#6a1b9a,stroke-width:2px 第1步:获取父本分子的预测分布:将原始的、未经修改的“父本分子”(parent molecule)输入到BCL-XpKa模型中,获取模型输出的pKa概率分布 Pparent。 第2步:对单个原子进行扰动:遍历分子中的每一个杂原子(非碳、氢原子),将其替换为一个保持价态和杂化状态正确的碳原子。 第3-4步:获取扰动分布并量化差异:将这个新的“扰动分子”输入到同一个BCL-XpKa模型中,获取其pKa概率分布 Pperturbed,并使用Kullback-Leibler(KL)散度来衡量 Pparent 和 Pperturbed 这两个概率分布的差异。 第5步:计算最终的ASA分数:原始的KL散度值需要经过一步经验性的去噪和放大,才能得到最终的ASA分数。其计算公式为: \(\text{ASA}(\text{atom } a) = e^{[S \cdot D_{KL}(P_{\text{perturbed}} || P_{\text{parent}})]} - 1\) 这个公式通过指数函数非线性地放大差异,使得影响显著的原子的分数远高于影响微弱的原子。 3. ASA的分析思路与应用 识别关键功能位点(Ionization Site Identification) 做法是直接找出分子中ASA分数最高的原子。这通常就是模型认为的、决定该分子pKa值的主要电离/质子化位点。论文在Novartis测试集上进行了验证,该测试集中的酸性分子平均有2.93个潜在电离亚结构,碱性分子平均有2.61个。ASA方法在识别最酸性原子时达到了灵敏度96.6%和特异性82.9%,在识别最碱性亚结构时也表现出色。这种方法的直接应用是在药物研发早期,可以快速、批量地为化合物库中的分子标注质子化状态,为后续的对接、MD模拟等步骤提供更准确的输入结构。 图3:用于分子分析的原子敏感性 (A) ASA协议示意图。“扰动pKa分布”和“父本pKa分布”指的是由BCL-XpKa输出的离散分布。 (B) 一个由BCL-XpKaAcid评分的酸的ASA分数示例。在这里,磺酰胺的氮原子被正确地选择为比其他潜在的酸性亚结构更具酸性。 (C) 一个由BCL-XpKaBase评分的碱的ASA分数示例。 (D) 一个碱,其中酰胺的氧原子在存在更具碱性的氮原子的情况下,主导了ASA分数。这种情况在61个含酰胺的碱性化合物中出现了4次。 (E) 用于BCL-XpKaBase分解的阳性(蓝色)和阴性(红色)对照亚结构的ASA分数。 (F) 通过添加一个胺基来调节吡啶氮的ASA分数,示例显示在x轴下方。ns = 不显著。*** = p < 0.001。 洞察模型的“化学知识”与学习机制(Probing Model Learning) 做法是比较同一官能团在不同化学环境下的ASA分数,或比较局部环境相同但整体化学性质迥异的基团。论文发现,邻近的吸电子基团(EWG)会显著降低吡啶氮的ASA分数(即降低其碱性贡献),而给电子基团(EDG)则会提高其分数(见图4F)。例如,在图2D的分子系列中,将哌啶(pKa 11.2,预测10.45)芳构化为吡啶(pKa 5.20,预测5.45),其碱性显著下降,模型准确地捕捉了这一趋势。这证明了模型虽然只学习了局部原子环境,但隐式地捕捉到了上下文依赖的化学规则。这种分析可以用来验证模型是否学到了正确的化学知识,而不是仅仅记住了某些表观特征,从而增加我们对模型预测的信任度。 图4:亚结构的原子敏感性分析 (A-D) 常见亚结构在作为分子的主要电离位点时与存在更主要电离位点时的ASA分数小提琴图。 (E) 常见亚结构在作为主要电离位点时的ASA分数小提琴图。 (F) 相邻的吸电子基团(EWG)和给电子基团(EDG)对吡啶氮ASA分数的影响,通过ASA分数的变化(ΔASA)来衡量。 (G) 分子对称性对ASA分数的“掩蔽效应”。 指导先导化合物的理性优化(Lead Optimization) 这是一个非常实用的应用场景,也是可以借鉴的完整工作流。论文以一个已知的KRAS降解剂PROTAC(P-1, PDB: 8QU8)为例,该分子与靶蛋白形成的复合物中,其连接臂上的一个叔胺与KRAS的Q62残基形成了关键的盐桥相互作用(图5A-B)。 流程: 问题定位:ASA分析显示,这个叔胺氮原子具有最高的ASA分数(12.1),确认了它就是导致PROTAC在生理条件下可能质子化的主要碱性位点(图5C)。而这种质子化状态不利于细胞膜通透性。 提出优化方案:针对这个高分子的叔胺,进行生物电子等排替换,例如将其替换为酰胺,设计出候选分子P-2(图5D)。 快速虚拟验证:BCL-XpKa模型预测P-2的pKa值显著降低至3.23,成功消除了碱性。 结构验证:通过对接模拟发现,新的P-2分子依然能够与KRAS的Q62残基形成一个关键的氢键,保持了必要的结合模式(图5F-H)。 结论:这个流程展示了如何利用ASA精确定位分子的物化性质缺陷来源,并指导进行高效、理性的化学修饰,从而在保持活性的前提下优化类药性。 图5:用于药物设计的原子敏感性分析 (A) 泛KRAS降解PROTAC P-1与VHL和KRAS形成的三元复合物的晶体结构(PDB: 8QU8)。 (B) PROTAC P-1,其pKa由BCL-XpKaBase计算为6.51。 (C) P-1连接臂中氮原子的ASA分数。 (D) 提出的P-1连接臂生物电子等排替换修饰及其由BCL-XpKa预测的pKa值。 (E) P-1和P-2对接到8QU8中VHL-KRAS蛋白-蛋白相互作用界面的三元复合物模型全局视图。 (F-H) 8QU8晶体结构和P-2酰胺修饰的结合位点视图,显示它们支持相似的PROTAC构象,并保留了与KRAS Q62的氢键。 三、 ASA的局限性与未来方向 作者在论文的讨论部分明确提到了当前框架的一些局限性,这对于我们借鉴和改进该方法至关重要。 仅限于原子级别:ASA旨在识别单个原子对预测的贡献,但不能直接输出“官能团级别”的贡献。一个原子的影响往往与它所在的整个官能团或药效团紧密相关,而ASA目前无法直接解耦这种集体效应。 对非直接电离原子的影响处理不完美:一些本身不电离但能显著影响pKa的原子(例如,通过强诱导效应或共振效应)偶尔会得到异常高的分数,从而干扰对真正电离位点的判断。例如,在一个含有酰胺的碱性分子中,有少数情况是酰胺的氧原子(本身不质子化)得到了最高的ASA分数,这可能是因为它被扰动后对分子整体电子云的改变过大,从而“掩盖”了真正的质子化位点。 扰动方式单一:将杂原子替换为碳是一种有意义但简化的扰动方式。对于某些复杂的化学环境,这种替换可能无法完全反映该原子在真实化学修饰中的作用。 未来方向: 指导数据增强:通过ASA识别出模型预测不佳或不确定的化学结构类型,可以指导性地扩充训练集,从而提升模型的性能和泛化能力。 整合到药物发现工作流:作者致力于将ASA整合到更大的药物发现工作流中,例如用于超大规模虚拟筛选(vHTS)的分子库预处理,以确保正确的分子质子化状态,提高筛选的命中率。 拓展到其他性质预测:论文提出,未来可以将ASA的思想应用于ADMET(吸收、分布、代谢、排泄、毒性)等更复杂性质的预测模型中,以理解和优化这些关键的药物属性。
Machine Learning & AI
· 2025-10-08
解构血脑屏障渗透性:一个可解释的多模态深度集成框架
解构血脑屏障渗透性:一个可解释的多模态深度集成框架 一、 论文整体概览 1. 论文基本信息 标题:Interpretable Multimodal Deep Ensemble Framework Dissecting Blood–brain Barrier Permeability with Molecular Features 中文译名:使用分子特征解构血脑屏障渗透性的可解释多模态深度集成框架 期刊:The Journal of Physical Chemistry Letters DOI: 10.1021/acs.jpclett.5c01077 发表年份:2025 Citation: J. Phys. Chem. Lett. 2025, 16, 5806-5819 2. 摘要、背景与科学问题 摘要翻译 血脑屏障渗透性(BBBP)预测在药物发现过程中扮演着关键角色,特别是对于靶向中枢神经系统(CNS)的化合物。尽管机器学习(ML)已显著推动了BBBP的预测,但目前仍迫切需要能够揭示调控BBB渗透性的物理化学原理的可解释性ML模型。在本研究中,我们提出了一个多模态ML框架,该框架整合了分子指纹(Morgan, MACCS, RDK)和图像特征以改进BBBP预测。分类任务(BBB可渗透 vs. 不可渗透)通过一个结合了多个基础分类器的堆叠集成模型来解决。在可比较的评估设置下,所提出的框架与近期的方法相比,展示了有竞争力的预测稳定性、泛化能力和特征可解释性。除了预测性能,我们的框架还结合了主成分分析(PCA)和沙普利加性解释(SHAP)分析,以突显对预测有贡献的关键指纹特征。回归任务(logBB值预测)则通过一个多输入深度学习框架来解决,该框架包含一个用于处理指纹的Transformer编码器,一个用于提取图像特征的卷积神经网络(CNN),以及一个用于增强特征交互的多头注意力融合机制。从多模态特征中提取的注意力图(Attention maps)揭示了分子表示内部的令牌(token)级关系。这项工作提供了一个可解释的框架,用于以增强的透明度和机理洞察力来建模BBBP,并为未来结合透明描述符和物理信息特征的研究奠定了基础。 背景与科学问题 血脑屏障(BBB)是保护中枢神经系统(CNS)的关键生理屏障,但它也成为CNS药物研发的巨大障碍。准确预测一个候选药物能否穿透BBB,是其成药性的决定性因素之一。传统的实验方法成本高昂且耗时,因此开发快速、准确且可靠的计算模型至关重要。 近年来,机器学习(ML)和深度学习(DL)模型在BBBP预测上取得了很高的准确率。然而,这些高性能的模型往往像一个“黑箱”,研究人员难以理解其做出特定预测的具体依据。这种可解释性的缺失不仅阻碍了我们对模型预测的信任,更重要的是,我们无法从模型学到的知识中提炼出清晰的、指导性的化学规则来辅助新药的理性设计。 因此,本文的核心科学问题是:如何在保证高预测精度的前提下,构建一个透明、可解释的BBBP预测框架,从而不仅“知其然”(预测结果),更能“知其所以然”(揭示分子结构与BBB渗透性之间的构效关系)? mindmap root(可解释性分析思路与实践) )特征集质量初评( ::icon(fa fa-flask) **PCA降维可视化** **核心思路**<br/>在建模之前快速评估特征集的质量与判别能力 **关键发现**<br/>MACCS指纹比Morgan指纹<br/>能更有效地分离BBB正负样本 **研究启发**<br/>这是筛选有效分子表示方法的一种重要且高效的前置步骤 )模型归因分析( ::icon(fa fa-search-plus) **SHAP值分析** **核心思路**<br/>定量计算每个“分子亚结构”对最终预测的贡献度 **高贡献度(高SHAP值)的关键亚结构** MACCS_43<br/>极性官能团(氢键供体/受体、磺酸盐) MACCS_39<br/>亚硫酸(酯) MACCS_37<br/>氨基酰胺(如脲结构) **化学原理验证** **极性表面积(PSA)原理**<br/>极性基团(如MACCS_43)增加PSA从而**降低**穿透脂质血脑屏障的能力<br/>(表现为高的负SHAP值) **卤化效应**<br/>MACCS_46(溴代基团)的负贡献可能源于分子量增加或代谢不稳定 **上下文依赖性** **结论**<br/>亚结构的最终效果受到<br/>分子整体拓扑和周围化学环境的共同调节 **具体表现**<br/>同一亚结构(如MACCS_38)在不同分子中可产生相反的SHAP贡献 **研究启发**<br/>为药物化学家提供**可操作的优化线索**<br/>指导基于构效关系的**理性药物设计** )深度模型内部机制探索( ::icon(fa fa-project-diagram) **注意力图可视化** **核心思路**<br/>揭示Transformer等深度模型在预测时“关注”的特征区域 **两种分析模式** **内部结构注意力**<br/>分析指纹序列内部各部分的重要性 **跨模态注意力**<br/>分析“指纹特征”与“图像空间区域”的关联 **关键发现**<br/>模型注意力会**从模糊逐渐聚焦**<br/>到化学上有意义的区域<br/>(例如`C=C`和`C-O-C`官能团) **研究启发**<br/>验证深度模型是否抓住了**正确的物理化学特征**<br/>为理解复杂模型的**内部工作机制**提供直观窗口 3. 模型框架总结 为解决上述问题,作者提出了一个多模态深度集成框架,其核心是融合不同来源的分子信息来提供更丰富的表征。 多模态特征输入:模型不依赖于单一的分子表示,而是同时整合了多种信息。其中分子指纹(Morgan、MACCS 和 RDK)和2D分子图像特征是两大核心输入模态。 指纹与模型的使用方式:论文中的框架分别针对分类和回归两个任务设计了不同的模型。在评估时,Morgan、MACCS和RDK这三种指纹是分开独立使用的,即用每一种指纹分别训练和评估模型,以比较不同分子表示方法的效果。它们并未融合成一个单一的巨大特征向量。 双任务模型架构:图2中展示的(a)和(b)是针对两个不同任务的两种独立模型。 分类模型(图2a,BBB+ vs. BBB-):该模型仅使用分子指纹作为输入。它采用一个堆叠集成模型(Stacking Ensemble Model)。该模型将多个基础分类器(如逻辑回归、随机森林、XGBoost等)的预测结果作为元特征(meta-features),再由一个最终的分类器进行综合决策,以提高模型的稳定性和泛化能力。 回归模型(图2b,logBB值预测):该模型采用了多模态输入,即同时使用分子指纹和2D分子图像。它是一个更复杂的多输入深度学习网络,使用Transformer处理序列化的指纹特征,用CNN处理图像特征,最后通过多头注意力机制(Multi-Head Attention)将这两种不同模态的特征进行深度融合。 图像特征的价值:论文明确提出,通过多模态融合来丰富分子表示是其核心策略之一。在回归模型中,作者专门设计了CNN和注意力模块来处理和融合图像特征。结论部分也强调,多样化分子模态的融合(即指纹+图像)结合透明的归因技术,能够提供更准确和有意义的预测。PCA分析(图8 e,f)显示,在与MACCS和RDK指纹融合后,特征空间的解释方差有所提升或保持高位,这表明图像特征确实为模型提供了有价值的互补信息,特别是在通过跨模态注意力分析揭示两种特征的关联时,其价值更为凸显。 二、 核心可解释性方法与发现 本文的亮点在于系统性地应用了多种前沿的可解释性技术来剖析其模型,从不同维度揭示了BBBP的分子层面的驱动因素。 1. 特征空间分析:PCA降维可视化 在进行复杂的模型解释之前,作者首先使用了主成分分析(PCA)这一经典的无监督降维方法,来直观地评估不同分子指纹对BBB+/BBB-两类分子的区分能力。 做法与发现: 将所有分子的Morgan指纹和MACCS指纹分别通过PCA降到二维空间进行可视化。 图6(a)显示,在使用Morgan指纹时,BBB+(红色)和BBB-(蓝色)两类分子的数据点严重重叠,难以区分,且前两个主成分仅能解释总方差的极小部分(PC1: 1.36%, PC2: 1.16%)。这表明Morgan指纹生成的特征向量虽然信息量大,但可能过于稀疏或其线性组合难以捕捉到类别间的清晰界限。 相比之下,图6(b)显示,在使用MACCS指纹时,两类分子的数据点形成了相对清晰可分的簇,且前两个主成分解释了更多的方差(PC1: 11.31%, PC2: 8.0%)。这说明MACCS指纹定义的166个预设化学亚结构,能够更有效地捕捉与BBB渗透性相关的结构差异。 应用与价值:PCA分析虽然简单,但它是在建模之前快速评估特征集质量和判别能力的有效手段。通过这种方法,作者在早期就得出结论:MACCS指纹在这种二元分类任务中比Morgan指纹更具信息量,这为后续选择MACCS作为主要特征进行SHAP分析提供了依据。 2. SHAP分析:量化分子指纹的贡献 SHAP(Shapley Additive Explanations)是一种源于合作博弈论的模型解释方法,它可以为单个样本的预测结果计算出每个输入特征的贡献值(SHAP值)。一个正的SHAP值表示该特征将预测推向正类(如BBB+),负值则推向负类(如BBB-)。 做法与发现: 作者对表现最好的分类模型(基于MACCS指纹的堆叠模型)进行了SHAP分析。图6(c)的蜂群图(Beeswarm plot)直观地展示了所有测试样本中,对模型影响最大的前几个MACCS指紋特征。 关键特征识别:分析发现,MACCS_43(通常代表富含氢键供体/受体和磺酸盐的极性官能团)、MACCS_39(O-S(=O)O)、MACCS_37(N-C(=O)N)和MACCS_38(N-C(-C)-N)等特征具有最高的平均SHAP值,表明它们对模型的预测有决定性影响。 化学意义的验证:图6(d)展示了包含这些关键亚结构的具体分子示例。例如,MACCS_43 在一个分子中贡献了+0.130的正SHAP值,而在另一个分子中贡献了-0.169的负SHAP值。这与化学直觉相符:极性基团通常会增加分子的极性表面积(PSA),从而降低其穿透富含脂质的血脑屏障的能力(对应负的SHAP值,预测为BBB-)。MACCS_46(代表溴代亚结构)的SHAP值为负,这可能是因为过度卤化会增加分子量或引入代谢不稳定性,从而整体上降低了渗透性。 上下文依赖性:SHAP分析还揭示了亚结构贡献的上下文依赖性。例如,MACCS_38(对称的二胺结构)在某些分子中贡献为正(+0.108),而在另一些分子中为负(-0.057),这表明其最终效果受到分子整体拓扑结构和周围化学环境的调节。 应用与价值:SHAP提供了一种强大的、定量的手段,可以将抽象的模型预测归因于具体的、化学家可以理解的分子亚结构。这使得模型的决策过程不再是“黑箱”,而是可以被验证和理解的。论文指出,这种由SHAP引导的分析为以结构功能关系为基础的CNS靶向药物设计提供了可操作的见解。虽然本文未直接展示用SHAP结果指导模型调优或新实验,但其揭示的关键特征无疑可以用于指导特征工程(例如,构建仅包含最重要特征的简化模型)或提出需要通过实验验证的化学假说(例如,系统性地修饰MACCS_43代表的基团来验证其对渗透性的影响)。 图6:特征分析与SHAP可解释性 (a) 使用Morgan特征的BBB+/BBB-样本的前两个主成分(PC)的PCA得分图。 (b) 使用MACCS特征的PCA得分图。 (c) 展示了使用MACCS指纹的模型中,样本级别的SHAP值分布的蜂群图。 (d) 包含关键亚结构的代表性分子的结构可视化,并标注了其对应的SHAP值(pos代表正贡献,neg代表负贡献)。 3. 注意力机制可视化:揭示模型内部焦点 对于用于logBB值预测的多输入深度学习模型,作者利用其核心组件——注意力机制(Attention Mechanism)——来探索模型在进行预测时,其“注意力”集中在哪些特征上。 做法与发现: 内部结构注意力(Intra-Structure Attention):图11(a)展示了Transformer编码器内部的注意力热图。这张图揭示了模型在处理一个分子的指纹序列时,不同指纹“令牌”(tokens)之间的相互依赖关系。 跨模态注意力(Cross-Modal Attention):图11(b)展示了从训练的第1个周期到第50个周期,分子指纹特征和CNN提取的图像特征之间的跨模态注意力图的演化。可以清晰地看到,随着训练的进行,模型逐渐学会将指纹中的特定信息(符号/化学特征)与图像中的特定空间区域(视觉特征)对应起来。模型的注意力从最初的模糊、分散状态,逐渐锐化并聚焦于化学上有意义的区域。 具体案例分析:以一个BBB+的分子divinylether为例,模型的高度注意力权重区域,无论是内部结构注意力还是跨模态注意力,都准确地对应于其分子结构中的关键官能团,如烯烃(C=C)和醚键(C-O-C)。 应用与价值:注意力可视化为理解深度学习模型(特别是基于Transformer的模型)的内部工作机制提供了一个直观的窗口。它能告诉我们模型在做决策时“正在看哪里”,从而验证模型是否抓住了正确的物理化学特征,而不是依赖于数据中的某些伪影或噪声。这种方法为以一种更具机理性的方式理解BBB渗透性提供了支持。 图11:注意力权重的可视化 (a) 代表指纹内部结构的注意力热图。 (B) 从第1个训练周期到第50个训练周期,结合图像特征的跨模态注意力热图的演化过程。 三、 本文的局限性与未来展望 作者在论文的结论部分坦诚地指出了当前工作的局限性,并对未来研究方向进行了展望。 需要更先进的特征选择技术:尽管当前框架表现良好,但未来可以引入更先进的特征选择方法(如LASSO、SHAP-RFE)来进一步优化输入特征,可能会提升模型性能和可解释性。 需要更广泛的外部验证:目前模型的验证主要基于B3DB数据集。为了证明其更广泛的适用性,未来需要在更多、更多样化的外部数据集上进行验证。 计算预测与实验验证的鸿沟:模型最终需要与真实的实验结果相结合。未来的工作需要整合实验分析,以弥合计算预测与药理学现实之间的差距。 展望:作者希望这个可解释的框架能够为药物发现早期阶段CNS活性化合物的设计和优先级排序做出贡献,并为未来融合更多透明化描述符和物理信息特征的研究铺平道路。
Machine Learning & AI
· 2025-10-08
<
>
Touch background to close