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
> QSAR application
A Bunch of Biophysics is Loading ...
QSAR application
QSAR算法排名在四大靶点上高度一致,但scaffold泛化差距因靶点而异
QSAR算法排名在四大靶点上高度一致,但scaffold泛化差距因靶点而异 本文信息 标题:系统性多靶点QSAR基准测试:机器学习算法、分子描述符与验证策略 作者:Salah A. Alshehade, Ghazi Al Jabal, Iqbal H. Jebril 发表期刊:Journal of Chemical Information and Modeling 发表时间:2026年(Received:2026年4月21日;Revised:2026年5月21日;Accepted:2026年5月26日) DOI:https://doi.org/10.1021/acs.jcim.6c01237 单位:Universiti Sultan Zainal Abidin(马来西亚)、MAHSA University(马来西亚)、Yarmouk Private University(叙利亚)、Al-Zaytoonah University of Jordan(约旦) 引用格式:Alshehade, S. A.; Al Jabal, G.; Jebril, I. H. Systematic Multi-Target QSAR Benchmarking: Machine Learning Algorithms, Molecular Descriptors, and Validation. J. Chem. Inf. Model. 2026. https://doi.org/10.1021/acs.jcim.6c01237 代码与数据:https://github.com/salahalsh/ML_QSARX (v1.0 tagged release);QSAR-X网页界面:https://insilicosigma.com/qsar-x/ 摘要 定量构效关系(QSAR)建模是计算药物发现的核心方法之一,但算法选择、描述符类型和验证策略对模型性能的影响尚未在多靶点、统一实验条件下得到系统性评估。本研究在四个治疗性靶点家族(EGFR激酶、DRD2 G蛋白偶联受体、BACE-1蛋白酶、hERG离子通道)的33751个化合物上,以完全一致的实验流程比较了10种机器学习算法和5种分子描述符。结果表明,scaffold划分导致的泛化差距在不同靶点间存在约2倍的变异($\Delta R^2$均值:0.084–0.171);适用域分析证实,超出Tanimoto化学域的化合物预测质量大幅下降($R^2$降幅0.31–0.51),其中hERG对结构远缘化合物的预测能力几乎丧失($R^2$:0.62 → 0.11);算法排名在四个蛋白家族间高度一致(Spearman $\rho$均值 = 0.92),树集成方法(随机森林 + ECFP4)在每个靶点上均优于基础图卷积网络(GCN,平均$R^2$亏损0.22)。 核心结论 算法排名跨靶点高度一致:随机森林、XGBoost、LightGBM的排序在四个靶点上几乎不变($\rho = 0.92$),说明最优算法选择具有可迁移性 scaffold泛化差距是数据属性,而非算法属性:hERG的泛化差距最大(0.171),DRD2和BACE-1最小(约0.085),且这一差距在不同算法间变化不大 ECFP4是最稳健的描述符选择:在全部四个靶点上表现最优或接近最优,优于ECFP6、MACCS和RDKit-2D 适用域报告不可或缺:域内$R^2$与域外$R^2$差距巨大,尤其对hERG而言,不报告AD等于掩盖模型的真实局限 基础GCN远不及精心构建的指纹基线:3层GCN(无边特征、无预训练)在所有靶点上均显著弱于RF + ECFP4 背景 药物发现是一项资源密集型工作,开发一种获批药物通常需要10–15年和超过26亿美元的投入。传统上,识别生物活性小分子主要依赖高通量筛选(HTS),但HTS成本高昂且受限于物理化合物库,而可药性分子空间的估计规模高达约 $10^{60}$ 个分子。定量构效关系(QSAR)建模通过将分子结构特征映射为生物活性预测值来应对这一挑战,其理论基础是Hansch-Fujita线性自由能框架。现代QSAR借助高维描述符和机器学习来捕获本质上非线性的结构-活性关系,已在虚拟筛选、先导化合物优化和ADMET预测中得到广泛应用。 在描述符层面:二维QSAR仅从分子图提取特征,无需三维坐标。扩展连通性指纹(ECFP)通过Morgan算法枚举原子周围的圆形化学环境,是当前最常用的分子表示之一;MACCS结构键以166位预定义子结构模式编码分子的宏观特征;RDKit-2D物理化学描述符则捕获分子量、logP、极性表面积等全局性质 在算法层面:随机森林(RF)因其鲁棒性强、对超参数不敏感、内部特征子采样天然适合高维稀疏指纹,已成为QSAR回归的事实标准。梯度提升方法(XGBoost、LightGBM、GBR)通过序列化误差校正可达到与RF相当的精度,而深度神经网络在表格型分子描述符数据上的提升并不一致 图神经网络的挑战:近年来,图神经网络(GNN)直接在分子图上端到端地学习任务特定的表示,理论上能捕获比手工描述符更丰富的结构信息,但在中等规模化学数据集(n < 10000)上是否能超越精心构建的指纹基线,仍缺乏严格对照 尽管算法和描述符选择的研究已有大量积累,但现有基准研究普遍存在几个共性问题: 不同研究使用的数据集、靶点、划分策略和评估指标差异巨大,跨研究的直接比较几乎不可行 大多数基准仅使用随机划分评估模型,而OECD QSAR验证原则强调的scaffold划分(Bemis-Murcko骨架)才能真正模拟对新骨架化合物的前瞻性预测能力 单靶点基准无法区分靶点特异性现象与普遍规律 本文通过在四个治疗性靶点家族(EGFR激酶、DRD2 GPCR、BACE-1蛋白酶、hERG离子通道)的33751个化合物上,以完全一致的实验流程系统比较10种算法和5种描述符,填补了这一空白。 已有基准研究的常见局限 本文的应对 单靶点评估,结论难以推广 四靶点跨家族(激酶/GPCR/蛋白酶/离子通道) 仅随机划分,$R^2$可能虚高 同时使用随机和scaffold划分,量化泛化差距 缺少适用域分析 Tanimoto距离AD分析 + 参数敏感性检验 描述符和算法比较不充分 10种算法 × 5种描述符的完整交叉对比 单次随机种子,结论不稳定 5个随机种子 + 1000次bootstrap置信区间 GNN对比缺少严格控制 3层GCN基线 vs RF + ECFP4,相同划分和数据 关键科学问题 算法选择是否具有跨靶点一致性? 在激酶、GPCR、蛋白酶和离子通道这四类截然不同的蛋白家族上,同一种算法是否始终表现最佳? 描述符层级的泛化性如何? ECFP4在文献中常被报告为最优描述符,但这一结论是否在多靶点、大样本量条件下依然成立? scaffold划分会暴露多大的泛化差距? 随机划分给出的乐观 $R^2$ 在多大程度上掩盖了模型在结构新颖化合物上的真实预测能力? 适用域(AD)能否量化预测可靠性? 当化合物超出训练集覆盖的化学空间时,模型性能下降多少?这一现象在不同靶点间是否一致? 创新点 四靶点统一流程基准测试:在涵盖激酶、GPCR、蛋白酶、离子通道的四个靶点上使用完全一致的实验流程,排除了方法学差异带来的混淆因素 多维度实验设计:原文方法部分列出9类系统实验,结果部分进一步以Experiment 10展开活性悬崖分析,覆盖算法比较、描述符比较、特征选择、划分策略、Y-scrambling验证、跨靶点一致性分析、超参数优化、GNN基线、适用域分析和活性悬崖分析 多随机种子验证与bootstrap置信区间:使用5个随机种子和1000次bootstrap重采样,确保结论不依赖于特定的数据划分 靶点依赖的泛化差距量化:在多靶点框架下系统报告scaffold gap的变异范围(约2倍),为QSAR模型的泛化能力评估提供了新的参考基准 研究内容 数据集与方法 研究通过ChEMBL REST API从ChEMBL 35数据库中提取了四个靶点的结合活性数据($\mathrm{IC}_{50}$或$K_i$),并使用自定义Python脚本(RDKit 2025.03.6)进行统一的数据清洗流程: 有效性过滤:去除缺少SMILES或活性值为非正的记录 去重:同一canonical SMILES对应多条测量值时取中位数 SMILES规范化:使用RDKit canonical SMILES标准化 盐去除:去除含片段分隔符的SMILES 分子量过滤:100–900 Da 活性转换:$\mathrm{IC}{50}$(nM)转为pActivity(= $-\log{10}(\mathrm{IC}_{50}/\mathrm{M})$) 最终获得33751个化合物,涵盖四个最重要的蛋白家族: 靶点 ChEMBL ID 蛋白家族 活性类型 原始记录 最终数量 pActivity范围 均值 ± SD EGFR CHEMBL203 激酶 $\mathrm{IC}_{50}$ 17652 10036 3.05–11.52 6.90 ± 1.35 DRD2 CHEMBL217 GPCR $K_i$ 13041 7558 3.07–11.52 6.77 ± 1.02 BACE-1 CHEMBL4822 蛋白酶 $\mathrm{IC}_{50}$ 14298 8080 3.00–12.00 6.68 ± 1.27 hERG CHEMBL240 离子通道 $\mathrm{IC}_{50}$ 12886 8077 3.00–9.85 5.43 ± 0.91 图1:四个靶点的pActivity分布直方图。EGFR(激酶,n=10036)和DRD2(GPCR,n=7558)的活性分布范围较宽,BACE-1(蛋白酶,n=8080)居中,hERG(离子通道,n=8077)的分布最窄(SD仅0.91),反映了hERG抑制剂的活性范围相对集中 描述符方面,研究计算了5种分子表示: ECFP4(半径=2,2048位):捕获中心原子周围2个化学键范围内的局部原子环境,是QSAR建模中最广泛使用的圆形指纹 ECFP6(半径=3,2048位):将捕获范围扩展至3个化学键,编码更大的分子子结构 MACCS结构键(166位):基于预定义的166个子结构模式的存在/不存在,是一种结构键指纹 RDKit 2D物理化学描述符(217维):通过RDKit的Descriptors.descList模块计算全部210个命名描述符(constitutional、topological、physicochemical),加上fragment-count描述符后保留217维。包括分子量、logP、拓扑极性表面积(TPSA)、氢键供体/受体数、可旋转键数等经典药物化学参数 组合描述符(RDKit-2D + ECFP4):将局部子结构信息与全局理化性质相结合 算法方面,研究比较了10种回归算法(详见下表),涵盖集成方法、核方法、实例学习、正则化线性模型和深度学习五大类。主体实验(Experiments 1-6, 8-10)使用下表默认超参数,Experiment 7专门对Top 4算法(RF、XGB、LGBM、SVR)进行了Bayesian超参数优化(Optuna TPE采样器,每靶点30次试验,5折交叉验证$R^2$目标),发现梯度提升方法获益最大(LGBM平均$\Delta R^2$ +0.040,XGB +0.034),而RF提升较小(+0.010)。 算法 类别 关键超参数 随机森林(RF) 集成 n_estimators=500, min_samples_split=5 XGBoost(XGB) 集成 n_estimators=300, max_depth=6, lr=0.1 LightGBM(LGBM) 集成 n_estimators=300, max_depth=6, lr=0.1 梯度提升(GBR) 集成 n_estimators=300, max_depth=5, lr=0.1 支持向量回归(SVR) 核方法 kernel=RBF, C=10, γ=scale K近邻(KNN) 实例学习 k=5, weights=distance Ridge 正则化线性 α=1.0 LASSO 正则化线性 α=0.01 弹性网络(Elastic Net) 正则化线性 α=0.01, l1_ratio=0.5 多层感知器(MLP) 深度学习 layers=(256,128), relu, adam, lr=0.001 graph TB direction TB subgraph S1["实验设计总览"] direction TB E1["1.算法比较:10种算法 × 5种描述符"] E2["2.描述符比较:RF × 5种描述符"] E3["3.特征选择:RDKit-2D逐步剪枝"] E4["4.划分策略:随机与Scaffold"] E5["5.Y-scrambling验证"] E6["6.跨靶点一致性:Spearman ρ"] E7["7.超参数优化:Optuna TPE"] E8["8.GNN基线:3层GCN"] E9["9.适用域:Tanimoto距离"] E10["10.活性悬崖:结构相似活性差异大"] end E1 --> E2 --> E3 --> E4 --> E5 E6 --> E7 --> E8 --> E9 --> E10 所有实验采用80:20的训练/测试划分。数据划分策略包括两种 orthogonal 方式: 随机划分:stratification-free随机采样,固定随机seed=42确保可复现性 Scaffold划分:使用RDKit的MurckoScaffold.GetScaffoldForMol()函数实现Bemis-Murcko scaffold划分 提取每个分子的Bemis-Murcko骨架(环系统+连接链) 无环化合物归为单一no-ring组 相同scaffold的分子聚成cluster 整个cluster只分配到训练集或测试集之一,确保scaffold不跨集 所有唯一scaffold cluster随机打乱(seed=42),贪心累积到测试集直到约20%总化合物数,剩余进训练集 单独出现的scaffold(singleton)作为独立cluster处理 评估指标包括决定系数$R^2$(coefficient of determination,衡量模型解释的方差比例)、RMSE和MAE,结果报告为5个随机种子(42、0、1、2、3)的均值 ± SD,并通过1000次bootstrap重采样计算95%置信区间。 实验1:算法比较——树集成方法的持续领先 以EGFR + ECFP4为例,10种算法在随机划分下的测试集$R^2$排名如下: 排名 算法 $R^2$(均值 ± SD) RMSE MAE 1 RF 0.726 ± 0.008 0.706 0.520 2 XGB 0.689 ± 0.011 0.752 0.574 3 SVR 0.692 ± 0.009 0.747 0.556 4 GBR 0.670 ± 0.010 0.774 0.594 5 LGBM 0.674 ± 0.010 0.769 0.587 6 Elastic Net 0.610 ± 0.017 0.841 0.648 7 LASSO 0.596 ± 0.014 0.857 0.664 8 KNN 0.604 ± 0.013 0.849 0.648 9 MLP 0.587 ± 0.017 0.866 0.637 10 Ridge 0.536 ± 0.028 0.917 0.695 RF以$R^2$ = 0.726排名第一,且bootstrap 95%置信区间证实RF对所有竞争算法的优势均具有统计显著性(CI不包含零)。前5名算法的$R^2$集中在0.67–0.73区间内,而线性模型(Ridge、LASSO)和MLP则明显落后。 图2:四个靶点上Top-5算法的测试集$R^2$比较(ECFP4指纹,随机划分)。每个子图顶部标注了跨种子排名的Spearman $\rho$值(均 > 0.92),表明算法排名在不同随机种子间高度稳定。RF在所有靶点上均排名第一或接近第一,其误差棒(SD)也是最小的 图3:EGFR上RF + ECFP4的预测值与观测值散点图。(a)训练集($R^2$ = 0.942):浅蓝色点沿红色虚线(完美预测线)紧密排列。(b)测试集($R^2$ = 0.720):数据点分散度增大,训练-测试差距为0.222,反映了RF在中等规模数据集上的典型过拟合程度。活性极端区域($\mathrm{pIC}_{50} < 5$和$> 10$)的偏差略有增大 图4:四个靶点上RF与Top竞争算法的$R^2$差异的bootstrap 95%置信区间(1000次重采样)。绿色柱表示RF优势具有统计显著性(CI不包含零),灰色柱标注“n.s.”表示差异不显著。RF在EGFR上对所有四个竞争算法均显著优于(绿色柱全为绿),在DRD2上对三个显著、对XGB不显著,在BACE-1上仅对SVR显著,在hERG上对所有算法均不显著 当将比较扩展到全部四个靶点时,这一模式高度一致:RF在EGFR和DRD2上明确排名第一;BACE-1上XGB以0.685略高于RF和LGBM(均为0.684),但bootstrap检验显示RF只显著优于SVR;hERG上RF和SVR同为0.578,所有Top算法差异均不显著。四个靶点的Top-5算法测试集$R^2$汇总如下(Table 4): 算法 EGFR DRD2 BACE-1 hERG RF 0.726 0.654 0.684 0.578 XGB 0.689 0.629 0.685 0.576 SVR 0.692 0.610 0.663 0.578 LGBM 0.674 0.611 0.684 0.557 GBR 0.670 0.612 0.673 0.559 图9:10种算法在4个靶点上的测试集$R^2$热力图。颜色从深红($R^2 \approx 0.3$)到深绿($R^2 \approx 0.8$),清晰展示了树集成方法(RF、XGB、LGBM、GBR)和SVR等强基线的优势。线性模型(Ridge)在所有靶点上排名末位 核心发现:树集成方法在QSAR回归任务上的统治地位不是偶然的。RF的内部特征bagging机制使其对高维稀疏特征(如2048位指纹)具有天然的鲁棒性,而梯度提升方法通过序列化学习也能有效处理这类数据。 实验2:描述符比较——ECFP4的稳健优势 使用RF作为固定算法,比较5种描述符在四个靶点上的表现: 描述符 EGFR DRD2 BACE-1 hERG ECFP4 0.726 0.654 0.684 0.578 ECFP6 0.720 0.641 0.678 0.577 MACCS 0.656 0.574 0.643 0.523 RDKit-2D 0.670 0.577 0.650 0.564 RDKit-2D+ECFP4 0.706 0.643 0.670 0.598 图5:5种描述符在4个靶点上的RF测试集$R^2$比较。ECFP4在EGFR、DRD2和BACE-1上均表现最佳;仅在hERG上,组合描述符(RDKit-2D + ECFP4)以0.020的$R^2$优势略胜。MACCS在四个靶点上整体最弱,尤其低于ECFP类指纹,说明166位结构键对精细局部化学环境的分辨力不足 描述符层级的稳定性值得注意:四个靶点上的描述符排名Spearman $\rho$均高于0.90,说明圆形指纹相对MACCS的优势并非特定于某个靶点,而是反映了ECFP在编码局部化学环境方面的固有信息优势。hERG的例外(组合描述符略优)提示,全局理化性质(如logP、极性表面积)对hERG阻断的预测提供了额外的互补信息。 实验3–4:特征选择与scaffold泛化差距 特征选择对RF的影响微乎其微。在EGFR上使用RDKit-2D描述符进行逐步特征剪枝:从无选择(217维)到方差过滤(183维)、相关性过滤( r ≤ 0.95,148维)、更严格的 r ≤ 0.75(112维),$R^2$仅从0.673下降至0.661。删去84.5%的特征仅损失0.012的$R^2$,这归功于RF的内部特征子采样机制。 图6:特征选择对EGFR上RF性能的影响。蓝色柱为测试集 $R^2$(左轴),橙色柱为保留的特征数(右轴)。五种特征选择配置下,$R^2$ 从0.673仅下降至0.661,删去84.5%的特征仅损失0.012 相比之下,scaffold划分揭示的泛化差距才是真正值得关注的问题。使用Bemis-Murcko scaffold划分后,$R^2$出现了靶点依赖的系统性下降: 靶点 RF XGB LGBM GBR 均值差距 EGFR 0.146 ± 0.061 0.145 ± 0.069 0.130 ± 0.057 0.139 ± 0.065 0.140 DRD2 0.102 ± 0.037 0.080 ± 0.042 0.080 ± 0.049 0.075 ± 0.038 0.084 BACE-1 0.097 ± 0.057 0.087 ± 0.044 0.080 ± 0.041 0.078 ± 0.040 0.085 hERG 0.179 ± 0.043 0.172 ± 0.065 0.158 ± 0.050 0.174 ± 0.047 0.171 图7:四个靶点的scaffold泛化差距(随机 $R^2$ 减去scaffold $R^2$),四种算法(蓝色=RF,绿色=XGB,橙色=LGBM,红色=GBR)在EGFR、DRD2、BACE-1、hERG上的表现。hERG的差距最大(均值0.171),DRD2和BACE-1最小(约0.085),呈现约2倍的变异范围。误差棒反映了5个随机种子的变异性 scaffold gap是数据集属性而非算法属性。四种算法在同一靶点上的差距高度一致,说明泛化困难源于数据本身的结构分布(如scaffold多样性、活性分布偏移),而非特定算法的过拟合倾向。 补充分析(Table S1)进一步揭示了scaffold划分导致的pActivity分布偏移:EGFR的测试集pActivity均值比训练集高0.249个log单位(KS检验p < 0.001),这是导致其较大scaffold gap的重要因素之一。 实验7:超参数优化——梯度提升获益最大 使用Optuna的TPE采样器进行贝叶斯超参数优化(每个靶点30次试验),三个规律浮现: 梯度提升方法获益最大:LGBM的$\Delta R^2$均值为+0.040,XGB为+0.034,而RF仅+0.010,SVR仅+0.007 优化后排名可能改变:XGB在优化后在四个靶点中的三个跃升至第一位,但这依赖于单次测试集划分,结论需谨慎对待 树方法收敛于狭窄的$R^2$区间:优化后,前5名算法的$R^2$差距不超过0.025,表明描述符信息量才是性能天花板 图10:默认超参数 vs Optuna优化后的测试集 $R^2$(30次Optuna试验,5折交叉验证)。四个子图分别为EGFR、DRD2、BACE-1、hERG,浅色柱为默认参数,深色柱为优化后参数。每个柱顶标注了$\Delta R^2$值。梯度提升方法(XGB、LGBM)的提升幅度显著高于RF和SVR 实验8:GNN基线——基础GCN不敌指纹 研究构建了一个3层图卷积网络(GCN),使用128维隐藏层、ReLU激活、全局平均池化,与RF + ECFP4进行对比: 靶点 划分方式 RF $R^2$ GCN $R^2$ $\Delta R^2$ EGFR 随机 0.720 0.513 −0.207 DRD2 随机 0.648 0.374 −0.274 BACE-1 随机 0.680 0.556 −0.124 hERG 随机 0.580 0.312 −0.268 EGFR Scaffold 0.474 0.275 −0.199 DRD2 Scaffold 0.541 0.321 −0.221 BACE-1 Scaffold 0.632 0.548 −0.084 hERG Scaffold 0.331 0.153 −0.178 图11:GCN vs RF + ECFP4在四个靶点上的 $R^2$ 对比。左图为随机划分,右图为scaffold划分。绿色柱为RF + ECFP4,紫色柱为3层GCN。红色数字标注了RF相对于GCN的优势(均为正值)。随机划分下DRD2的差距最大($\Delta R^2$ = −0.274),scaffold划分下DRD2的差距也最大($\Delta R^2$ = −0.221) GCN在随机划分下平均亏损0.218 $R^2$,在scaffold划分下平均亏损0.170 $R^2$。随机划分下差距最大的靶点是DRD2(亏损0.274),hERG也接近这一水平(亏损0.268)。研究使用的GCN不含边特征和预训练,属于最基础的图架构,但即便如此,这一结果也提醒我们:在中等规模的化学数据集(约7500–10000个化合物)上,GNN的数据效率远不及精心构建的分子指纹。 实验9:适用域——预测可靠性的量化 基于Tanimoto距离的适用域分析(k=5近邻,Z=1.5阈值)揭示了模型预测的可靠性边界: 靶点 AD阈值 域内占比 $R^2$(域内) $R^2$(域外) $R^2$降幅 EGFR 0.475 92.7% 0.738 0.369 −0.370 DRD2 0.468 92.0% 0.672 0.322 −0.350 BACE-1 0.447 92.3% 0.680 0.366 −0.314 hERG 0.590 87.4% 0.620 0.108 −0.512 图12:Williams图——四个靶点的标准化残差(y轴)与Tanimoto距离(x轴)关系。蓝色圆点为域内化合物(Inside AD),红色叉号为域外化合物(Outside AD),红色虚线为AD边界。各靶点的AD阈值分别为:EGFR 0.475、DRD2 0.468、BACE-1 0.447、hERG 0.590。域外化合物的残差离散度显著高于域内 图13:域内与域外预测质量对比。灰色柱为全部测试集 $R^2$,绿色柱为域内 $R^2$,红色柱为域外 $R^2$。hERG的域外 $R^2$ 仅0.107,几乎丧失预测能力 hERG的适用域问题最为严峻:约13%的测试化合物位于AD之外,这些化合物的$R^2$从0.62骤降至0.11,意味着模型对结构远缘的hERG抑制剂几乎丧失了预测能力。这一发现对hERG安全性预测的实际应用提出了重要警示。 实验10:活性悬崖——描述符的分辨力差异 活性悬崖(Activity Cliff)是指结构高度相似但活性差异巨大的化合物对,是QSAR建模的“天敌”。研究定义cliff对的条件为ECFP4 Tanimoto $\ge 0.6$且$ \Delta \mathrm{pActivity} \ge 2.0$: BACE-1的cliff密度最高(9531对,39.6%的化合物;测试集中38.2%的化合物属于cliff),其次是EGFR(8547对,37.0%)。EGFR上ECFP4的cliff化合物$R^2$为0.632,非cliff化合物为0.773 hERG的cliff密度最低(10.8%),且cliff与非cliff的$R^2$差距也最小($\Delta = 0.100$),反映了hERG数据集较低的活性异质性 MACCS的cliff预测能力最差:四个靶点的MACCS cliff $R^2$均低于ECFP4和ECFP6,尤其在DRD2上仅为0.316,因为166位结构键指纹无法区分仅差几个原子的cliff对 图15:描述符类型对活性悬崖化合物预测的影响(RF,随机划分,seed=42)。四个面板分别为EGFR、DRD2、BACE-1、hERG。每个面板中,蓝色柱为非cliff化合物的 $R^2$,灰色柱为全部测试化合物的 $R^2$,红色柱为cliff化合物的 $R^2$。Cliff对定义:ECFP4 Tanimoto $\ge 0.6$且$ \Delta \mathrm{pActivity} \ge 2.0$。各面板上方标注了cliff和非cliff化合物的数量。ECFP4在EGFR、DRD2和BACE-1上给出最高cliff $R^2$,hERG上则是RDKit-2D + ECFP4略高 残差分析(图14,见附录)证实了RF + ECFP4模型在EGFR上无系统性偏差,MAE = 0.520 $\mathrm{pIC}_{50}$单位,残差分布近似对称。 关键结论与批判性总结 本文在多靶点、统一流程下对QSAR建模的关键决策因素进行了系统性评估,主要贡献包括: 为QSAR建模提供了可操作的决策指南:RF + ECFP4作为默认配置在绝大多数场景下是合理的选择;超参数优化优先应用于梯度提升方法 量化了scaffold泛化差距的靶点依赖性:约2倍的变异范围(0.084–0.171)提醒研究者,不同靶点的结构外推难度可能截然不同 强调了适用域报告的必要性:域内外$R^2$差距高达0.31–0.51,不报告AD等于对使用者隐瞒了模型的关键局限 存在的局限性: 特征选择仅在EGFR上进行了详细分析,结论的跨靶点普适性有待验证 超参数优化只使用30次Optuna试验和单一测试划分,优化后XGB排名上升这一现象仍需要多随机种子确认 ECFP指纹折叠为2048位向量,哈希折叠不可避免地带来一定bit collision风险 GNN基线仅测试了最简单的GCN架构,未包含边特征、attention、预训练或Chemprop、AttentiveFP、Uni-Mol等更强模型 ChEMBL assay异质性仍然存在,同一化合物多条记录取中位数只能部分缓解实验噪声 本文限于二维描述符和四类靶点,三维描述符、interaction fingerprint、physics-informed表示以及更多靶点类别仍需要进一步验证
Machine Learning & AI
· 2026-06-22
附录:QSAR基准测试的技术细节与补充分析
附录:QSAR基准测试的技术细节与补充分析 本附录补充主文档中未展开的技术细节,包括Y-scrambling验证、bootstrap置信区间分析、scaffold划分的详细统计量、适用域敏感性分析、EGFR残差分析,以及活性悬崖的完整描述符对比。 Y-scrambling验证 Y-scrambling(响应变量随机打乱)是验证QSAR模型是否捕获了真实的构效关系信号、而非偶然相关性的标准方法。本研究在LightGBM + ECFP4上进行了30次打乱迭代: 靶点 原始CV $R^2$ 打乱均值 打乱SD 差距 EGFR 0.654 −0.069 0.007 0.723 DRD2 0.603 −0.084 0.010 0.687 BACE-1 0.670 −0.079 0.010 0.749 hERG 0.516 −0.086 0.010 0.602 图8:四个靶点的Y-scrambling验证。红色直方图为30次打乱后的CV $R^2$分布(均值接近零或为负值),红色虚线为原始模型的CV $R^2$。原始模型与打乱分布之间的巨大差距(0.60–0.75 $R^2$)证实了模型捕获的是真实的结构-活性关系信号 所有四个靶点的打乱$R^2$均值均在零附近或为负值(−0.069至−0.086),与原始模型$R^2$的差距均超过0.60。这一结果排除了算法偏差或数据伪影导致虚假高性能的可能性。 Bootstrap置信区间与算法差异的统计检验 为了评估RF对其他算法的优势是否具有统计显著性,研究计算了1000次bootstrap重采样的95%置信区间。以EGFR为例(ECFP4,随机划分): 对比 $\Delta R^2$(RF − 对手) 95% CI 是否显著 RF vs XGB +0.039 [+0.028, +0.050] 是 RF vs LGBM +0.051 [+0.039, +0.062] 是 RF vs SVR +0.031 [+0.012, +0.049] 是 RF vs GBR +0.061 [+0.048, +0.076] 是 图3:EGFR上RF + ECFP4的预测值与观测值散点图。训练集$R^2$ = 0.942(a),测试集$R^2$ = 0.720(b),训练-测试差距为0.222,反映了RF在中等规模数据集上的典型过拟合程度。测试集中数据点沿完美预测线的分散程度在活性范围两端略有增大,这与训练数据密度效应一致 图4:四个靶点上RF与Top竞争算法的$R^2$差异的bootstrap 95%置信区间。绿色柱表示RF优势具有统计显著性(CI不包含零),灰色柱表示差异不显著。RF在EGFR和DRD2上的优势较稳定,但在BACE-1和hERG上的多数小幅差异并不显著 统计显著性具有明显的靶点依赖性:RF在EGFR上相对四个竞争算法均显著占优,在DRD2上相对LGBM、SVR和GBR显著但相对XGB不显著;在BACE-1上仅相对SVR显著;在hERG上所有对比均未达到显著。这说明RF是稳健强基线,但不能把所有小幅$R^2$差异都解释为真实优势。 Scaffold划分的详细统计量 Bemis-Murcko scaffold划分(seed = 42,目标80:20比例)的详细统计信息: 靶点 总化合物 唯一scaffold 孤立scaffold(%) 最大cluster 训练集 测试集 测试占比 pActivity偏移(KS检验) EGFR 10036 3700 2490(67.3%) 547 7780 2256 22.5% +0.249(p < 0.001) DRD2 7558 3549 2501(70.5%) 69 6036 1522 20.1% +0.008(p = 0.256) BACE-1 8080 3156 2089(66.2%) 146 6463 1617 20.0% +0.168(p < 0.001) hERG 8077 4131 2977(72.1%) 103 6460 1617 20.0% +0.086(p < 0.001) 几个关键观察: EGFR的scaffold cluster最大(547个化合物),且scaffold多样性较低(3700个唯一scaffold对应10036个化合物),这与其激酶抑制剂的高度保守骨架结构一致 DRD2的pActivity偏移最小(+0.008,KS p = 0.256),说明训练集和测试集在活性分布上高度一致,其较小的scaffold gap(0.084)不能归因于分布偏移 EGFR的pActivity偏移最大(+0.249,KS p < 0.001),测试集的活性均值高于训练集,这意味着scaffold划分将更多高活性化合物分入了测试集,这是其较大scaffold gap的重要贡献因素之一 hERG的孤立scaffold比例最高(72.1%),反映了hERG抑制剂化学空间的高度碎片化 适用域敏感性分析 AD阈值由两个参数控制:k近邻数$k$和Z-score阈值$Z$。默认设置为$k=5$、$Z=1.5$。以下展示了参数变化对覆盖率和$R^2$降幅的影响(RF + ECFP4,随机划分): 靶点 k Z=1.0(覆盖率 / $R^2$降幅) Z=1.5(覆盖率 / $R^2$降幅) Z=2.0(覆盖率 / $R^2$降幅) EGFR 3 88.7% / 0.298 93.1% / 0.454 95.6% / 0.697 EGFR 5 88.2% / 0.285 92.7% / 0.370 95.4% / 0.729 EGFR 10 87.0% / 0.228 92.0% / 0.348 95.2% / 0.630 DRD2 5 85.7% / 0.209 92.0% / 0.350 96.2% / 0.457 BACE-1 5 87.1% / 0.171 92.3% / 0.314 95.5% / 0.461 hERG 5 82.1% / 0.465 87.4% / 0.512 93.7% / 0.528 hERG对AD参数的变化最为敏感:即使在最宽松的设定(k=5, Z=1.0)下,hERG的$R^2$降幅仍高达0.465。而在最严格的设定(k=5, Z=2.0)下,hERG的域外$R^2$仅比0高出约0.05。这一敏感性分析的定性结论(hERG > EGFR > DRD2 ≈ BACE-1)在所有参数组合下保持一致,增强了结论的可靠性。 残差分析 图14:EGFR上RF + ECFP4的残差分析。(a)残差 vs 预测值散点图:x轴为预测$\mathrm{pIC}_{50}$(4–10),y轴为残差(观测 − 预测),红色虚线为零残差参考线。数据点(浅蓝色)围绕零线随机分布,无系统性偏差。(b)残差分布直方图(MAE = 0.520):蓝色柱为频率分布,红色虚线为均值残差位置,分布近似对称 EGFR上RF + ECFP4(随机划分)的残差分析显示: 无系统性偏差:残差在pActivity预测值4–10的范围内围绕零线随机分布,无明显的趋势或模式 MAE = 0.520 $\mathrm{pIC}_{50}$单位,即平均预测误差约半个log单位 残差分布近似对称:直方图呈近似正态分布,峰值在零附近 活性极端区域的方差略有膨胀:在pActivity < 4.5和 > 10.5的区域,残差离散度增大,这与训练数据密度在这些区域的稀疏性一致 活性悬崖的完整描述符对比 完整的活性悬崖分析结果(RF,随机划分,seed = 42),cliff对定义为ECFP4 Tanimoto $\ge 0.6$且$ \Delta \mathrm{pActivity} \ge 2.0$: 靶点 描述符 cliff测试化合物 非cliff测试化合物 $R^2_{all}$ $R^2_{cliff}$ $R^2_{non-cliff}$ EGFR ECFP4 765 1243 0.720 0.632 0.773 EGFR ECFP6 765 1243 0.716 0.631 0.765 EGFR MACCS 765 1243 0.646 0.543 0.706 EGFR RDKit-2D 765 1243 0.671 0.575 0.727 EGFR RDKit-2D+ECFP4 765 1243 0.706 0.613 0.761 DRD2 ECFP4 219 1293 0.641 0.476 0.682 DRD2 ECFP6 219 1293 0.634 0.449 0.681 DRD2 MACCS 219 1293 0.541 0.316 0.597 DRD2 RDKit-2D 219 1293 0.556 0.375 0.599 DRD2 RDKit-2D+ECFP4 219 1293 0.628 0.437 0.677 BACE-1 ECFP4 618 998 0.682 0.517 0.767 BACE-1 ECFP6 618 998 0.679 0.517 0.761 BACE-1 MACCS 618 998 0.638 0.457 0.731 BACE-1 RDKit-2D 618 998 0.650 0.465 0.746 BACE-1 RDKit-2D+ECFP4 618 998 0.671 0.495 0.763 hERG ECFP4 174 1442 0.584 0.481 0.581 hERG ECFP6 174 1442 0.575 0.458 0.577 hERG MACCS 174 1442 0.550 0.426 0.551 hERG RDKit-2D 174 1442 0.568 0.478 0.557 hERG RDKit-2D+ECFP4 174 1442 0.595 0.497 0.591 ECFP4在EGFR和DRD2两个靶点的cliff化合物预测上表现最佳,在BACE-1上与ECFP6并列最高(均为0.517);hERG上则是RDKit-2D + ECFP4组合描述符最高(0.497),略高于ECFP4(0.481)。MACCS的cliff预测$R^2$始终最低,尤其是对DRD2($R^2_{cliff}$仅0.316),进一步印证了预定义结构键指纹在精细结构区分上的局限性。 hERG的cliff密度最低(2529对,占全部化合物10.2%;测试集中174个cliff化合物,占10.8%),且cliff与非cliff化合物的$R^2$差距也最小(ECFP4:$\Delta R^2$ = 0.100),这与其较窄的pActivity分布(SD = 0.91)一致——活性范围有限意味着“cliff”($ \Delta \mathrm{pActivity} \ge 2.0$)的绝对数量较少。
Machine Learning & AI
· 2026-06-22
整合qHTS与QSAR:筛选hERG风险较低的GPCR先导化合物
整合qHTS与QSAR:筛选hERG风险较低的GPCR先导化合物 本文信息 标题:整合qHTS和QSAR模型以识别安全的GPCR靶向化合物:关注hERG依赖性心脏毒性 作者:Xi Luo, Jinghua Zhao, Srilatha Sakamuru, Menghang Xia, Tuan Xu, Jameson Travers, Carleen Klumpp-Thomas, Hu Zhu, Matthew D. Hall, Stephen S. Ferguson, David M. Reif, Ruili Huang 发表时间: 2026年2月17日 单位: 美国国家推进转化科学中心(NCATS)、北卡罗来纳大学等(美国) 引用格式: Luo, X., Zhao, J., Sakamuru, S., Xia, M., Xu, T., Travers, J., Klumpp-Thomas, C., Zhu, H., Hall, M. D., Ferguson, S. S., Reif, D. M., & Huang, R. (2026). Integrating qHTS and QSAR Models to Identify Safe GPCR-Targeted Compounds: A Focus on hERG-Dependent Cardiotoxicity. Journal of Chemical Information and Modeling, 66(7), 2474–2487. https://doi.org/10.1021/acs.jcim.5c02291 相关工具:ChemoTyper(ToxPrint chemotypes)https://github.com/mn-am/chemotyper 摘要 G蛋白偶联受体是七跨膜受体家族,通过G蛋白介导细胞外信号转导,在多种生理和神经过程中发挥关键作用。ADRB2、CHRM1、DRD2和HTR2A等重要GPCR靶点,与哮喘、精神分裂症等疾病的治疗密切相关。然而,许多靶向GPCR的药物会抑制hERG钾离子通道,导致QT间期延长,也就是心电图上反映心室去极到复极全过程的时间变长,并增加心律失常风险。本研究整合定量高通量筛选和基于机器学习的定量结构活性关系模型,采用不同的数据处理顺序预测hERG风险较低的选择性GPCR靶向化合物。模型在Tox21 10K化合物库上训练,经LOPAC数据集,即Library of Pharmacologically Active Compounds,外部验证,随后用于虚拟筛选约36万个多样化化合物,并对预测排名靠前的化合物进行实验验证,发现了多个hERG风险较低的新型GPCR调节剂。 核心结论 hERG毒性普遍存在:在GPCR活性化合物中,尤其是拮抗剂模式中,hERG抑制剂的占比接近或超过50%的GPCR活性拮抗剂,强调在GPCR药物开发中监测hERG抑制的重要性 双模型策略有效:Model 1和Model 2都能给出稳定预测,最优模型的AUC-ROC可达AUC-ROC值0.84以上 共识模型成功筛选:使用四种ML算法(RF、SVM、NB、XGB)的共识策略筛选1408个CHRM1预测活性化合物,实验验证显示激动剂PPV达阳性预测值0.72,拮抗剂PPV达阳性预测值0.91,hERG模型的NPV为阴性预测值81.6% 发现新型先导化合物:鉴定出多个具有微摩尔级活性的CHRM1激动剂和拮抗剂,且hERG抑制较弱,说明这套流程适合用于早期候选物优先级排序 背景 G蛋白偶联受体是最大的细胞表面受体家族,跨越细胞膜七次,通过细胞外环与配体结合,通过细胞内环与G蛋白相互作用。GPCR在各种生理和神经过程中至关重要,是哮喘、阿尔茨海默病、帕金森病、精神分裂症等多种疾病的治疗靶点。例如,β2肾上腺素受体激动剂如沙丁胺醇用于治疗哮喘,毒蕈碱乙酰胆碱受体M1激动剂如占诺美林可改善阿尔茨海默病的认知功能,多巴胺D2受体激动剂如普拉克索用于帕金森病,5-羟色胺受体2A拮抗剂如氯氮平用于精神分裂症。 然而,许多靶向GPCR的药物与心脏毒性副作用相关,这主要归因于它们对hERG(human Ether-à-go-go-Related Gene)钾离子通道的抑制作用。hERG编码Kv11.1,是延迟整流钾通道快速组分的α亚基,对心脏复极化至关重要。抑制hERG通道会导致QT间期延长。这里的QT间期,指的是心电图中从Q波起点到T波终点的一段时间,可粗略理解为心室完成一次电活动所需的时间。这个时间一旦拉长,就会增加尖端扭转性室性心动过速等严重心律失常的风险,可能进展为室颤和猝死。因此,hERG抑制是药物淘汰和市场撤市的主要原因,FDA要求几乎所有新的低分子量药物都必须进行“全面QT”研究以评估其对QT间期延长的影响。在药物开发早期识别hERG抑制对于预防心脏毒性、提高药物安全性、确保监管合规和优化药物开发过程至关重要。 定量高通量筛选是一种强大的工具,可用于识别各种分子靶点的新型先导化合物。Tox21计划应用qHTS测试了约10000个药物和环境化学物质(Tox21 10K化合物库),涵盖约80个体外实验,包括核受体、应激反应通路、GPCR以及其他毒性相关靶点。重要的是,扩展的Tox21实验组合还包括专门的hERG通道抑制实验,提供了关键心脏毒性终点的直接测量。Tox21实验数据已用于构建毒性预测模型以及识别疾病靶点的新型先导化合物。 基于机器学习的定量结构活性关系模型是传统湿实验室实验的实用且有效的替代方案,已被用于虚拟筛选大型化学库,以识别GPCR激动剂、拮抗剂以及hERG抑制剂。qHTS实验数据为开发ML模型提供了稳健的数据集,用于预测小分子对不同靶点,如GPCR与hERG的活性和选择性。先前研究已经报道,ML模型可以成功识别具有GPCR活性和hERG抑制活性的分子。然而,设计用于识别GPCR活性化合物的机器学习模型也可能同时选出抑制hERG的候选物。因此,需要在药物发现早期优先考虑兼具GPCR活性和较低hERG风险的虚拟筛选方法。 关键科学问题 GPCR药物的心脏毒性风险:许多靶向GPCR的药物会抑制hERG通道,导致QT间期延长和心律失常,如何在药物开发早期有效识别和排除hERG抑制剂? 选择性预测的挑战:如何构建能够同时预测GPCR活性和hERG抑制的机器学习模型,以筛选出具有选择性的安全先导化合物? 数据不平衡问题:在GPCR活性化合物中,hERG抑制剂的比例很高(尤其是拮抗剂),如何处理这种数据不平衡并训练稳健的分类模型? 模型泛化能力:如何确保模型在化学结构多样的化合物库中保持良好的预测性能,并成功应用于外部验证和大规模虚拟筛选? 创新点 双模型策略:提出两种不同的建模策略,Model 1分别为8个GPCR靶点和hERG构建独立模型,Model 2在构建GPCR模型前排除hERG抑制剂,系统比较了两种策略的性能 整合qHTS与QSAR:利用Tox21 10K化合物库的qHTS数据构建ML模型,结合ECFP4指纹和多种ML算法,实现了从高通量筛选数据到虚拟筛选的有效转化 共识模型筛选:采用四种ML算法(RF、SVM、NB、XGB)的共识策略筛选约36万个化合物,并通过严格的hERG排除阈值0.3(预测概率≥预测概率阈值0.3)降低心脏毒性风险 实验验证成功:对模型预测的CHRM1活性化合物进行实验验证,发现了多个具有微摩尔级活性且无明显hERG抑制的新型先导化合物,验证了模型的实用性 研究内容 本研究整合定量高通量筛选和机器学习QSAR模型,旨在开发能够预测选择性GPCR靶向化合物,即hERG风险较低候选物的计算方法。研究针对四个重要的GPCR靶点,即ADRB2、CHRM1、DRD2和HTR2A的激动剂和拮抗剂模式,采用两种不同的数据建模流程,即Model 1和Model 2构建分类模型,通过Tox21 10K化合物库的qHTS数据训练,LOPAC数据集外部验证,最终应用于NCATS内部约36万个化合物的虚拟筛选,并对预测排名靠前的化合物进行实验验证。 方法详述 数据来源 Tox21 10K化合物库包含8599个独特化合物,其中约3000个为获批药物。研究通过qHTS获得四个GPCR靶点,即ADRB2、CHRM1、DRD2和HTR2A的激动剂与拮抗剂活性数据,以及hERG通道抑制数据。每个化合物都在15个浓度下进行三重复测试。 数据处理流程 曲线分级:根据浓度-响应曲线观察到的形状分配类别(1.1-1.4和2.1-2.4为活性,3为活性,4为非活性) 曲线秩次:转换为-9到9之间的整数,秩次越高表示曲线质量、效力和有效性越高。抑制剂分配负秩次,激活剂分配正秩次 活性判定:基于平均曲线秩次和三次重复实验的重现性,将化合物分配为“活性激动剂/拮抗剂”、“非结论性激动剂/拮抗剂”、“非结论性”或“非活性” 图1:模型构建和外部验证的数据集与框架 图1A:hERG抑制剂(橙色段)在八个数据集的活性GPCR化合物中的分布(包含橙色和蓝色段的柱子),包括ADRB2、CHRM1、DRD2和HTR2A的激动剂和拮抗剂 图1B:GPCR实验数据中活性和非活性化合物的分布 图1C:两种ML模型使用的数据集中活性和非活性化合物的分布,Model 1分别为8个GPCR靶点和hERG构建独立模型,Model 2从GPCR活性化合物中排除hERG抑制剂 图1D:外部验证数据集(LOPAC)中活性和非活性化合物的分布 图1E:虚拟筛选约36万个多样化化合物并对选定的预测进行实验验证的流程 这张图把整篇文章的逻辑压缩得很清楚。图1A先说明问题本身,即活性GPCR化合物里混有大量hERG抑制剂;图1C再展示两种建模流程的差别;图1D和图1E则对应外部验证与大规模虚拟筛选,基本就是全文的方法主线。 双模型建模策略 本研究采用两种不同的建模流程来预测选择性GPCR靶向化合物: graph TB subgraph S2["Model 2:预先排除策略"] direction TB B1["Tox21 10K数据集"] B2["识别hERG抑制剂<br/>并从GPCR活性中排除"] B3["8个GPCR数据集<br/>已排除hERG抑制剂"] B4["训练8个GPCR模型<br/>只包含非hERG抑制剂的GPCR活性化合物"] B5["预测GPCR活性<br/>默认低hERG风险"] B1 --> B2 B2 --> B3 B3 --> B4 B4 --> B5 end subgraph S1["Model 1:分别建模策略"] direction TB A1["Tox21 10K数据集"] A2["8个GPCR数据集<br/>ADRB2/CHRM1/DRD2/HTR2A<br/>激动剂+拮抗剂"] A3["hERG抑制剂数据集"] A4["独立训练9个模型<br/>8个GPCR模型 + 1个hERG模型"] A5["独立预测<br/>GPCR活性 + hERG抑制"] A1 --> A2 A1 --> A3 A2 --> A4 A3 --> A4 A4 --> A5 end Model 1采用分别建模策略,为8个GPCR靶点和hERG构建独立的分类模型,优点是灵活性高,可根据实际需求调整GPCR活性和hERG毒性的权重。Model 2采用预先排除策略,在训练GPCR模型前先排除hERG抑制剂,直接训练选择性模型,优点是简化后续筛选流程。通过对比两种策略,可以系统评估先识别活性、再剔除hERG风险与直接训练选择性模型的优劣。 分子描述符 ECFP4(Extended Connectivity Fingerprints 4)为1024位指纹,编码局部原子环境,如原子类型、芳香性、环成员、杂原子和键序,用来捕获常见亚结构特征。 机器学习算法 算法 作用特点 朴素贝叶斯 概率分类器,假设特征之间相互独立 随机森林 集成学习方法,通过多棵决策树投票得到结果 支持向量机 通过寻找最优超平面拉开不同类别间隔 XGBoost 梯度提升树方法,迭代优化分类误差 模型评估 项目 设置 交叉验证 5折分层交叉验证,重复10次 性能指标 AUC-ROC、平衡准确率、马修斯相关系数 类别平衡 在训练集上使用随机欠采样 共识策略 使用四种经过验证的机器学习分类器,即RF、SVM、NB和XGB,在Tox21 10K化合物库上训练并经LOPAC数据集外部验证的模型,对NCATS内部约36万个化学多样性化合物进行虚拟筛选。如果四个模型独立给出的活性概率都高于各自阈值,化合物才会被判定为GPCR活性。 hERG排除 为最大限度降低心脏毒性风险,研究统一使用hERG排除阈值0.3:凡是预测hERG抑制概率大于等于阈值0.3的化合物都会被排除。由于资源限制,最终每个GPCR靶点只保留约2000个候选,优先进入实验的是预测GPCR活性更高、预测hERG风险更低的那一批。 实验验证 基于四种ML模型的预测概率,研究选择模型预测的CHRM1活性化合物进行实验验证。总计测试1408个化合物,其中包括382个预测激动剂和1037个预测拮抗剂,另有12个化合物同时被预测为激动剂与拮抗剂。这些样品随后在CHRM1激动剂模式、CHRM1拮抗剂模式和hERG抑制实验中接受测试。 结果与分析 hERG毒性在GPCR药物中的普遍性 图1A揭示了hERG抑制剂在GPCR活性化合物中的广泛分布。例如,在45个ADRB2活性激动剂中,有13个化合物是hERG抑制剂。在其他GPCR活性化合物中也发现了大量的hERG抑制剂,尤其是在拮抗剂模式实验中,接近或超过50%的GPCR活性拮抗剂也抑制hERG。这种高比例的hERG毒性表明,单纯筛选GPCR活性化合物不足以确保药物安全性,必须同时评估hERG抑制风险。 模型训练性能评估 图2:Model 1(左)和Model 2(右)的性能 使用四种ML算法(NB、RF、SVM和XGB)开发的模型通过受试者工作特征曲线下面积(AUC-ROC)、平衡准确率和马修斯相关系数进行评估 指标报告为10次5折分层交叉验证中各折的平均值±标准差 在每一折中,数据集分为训练和测试子集,对训练数据应用随机欠采样以处理类别不平衡,并通过评估预测概率与测试集对比来计算AUC-ROC、BAC和MCC指标 图2的重点不是某一个单独柱子有多高,而是两个关键观察。 第一,不同算法之间确实有差异,但多数任务都能维持在可用区间,说明数据本身足以支撑分类建模。 第二,Model 2在大多数GPCR任务上的AUC-ROC略高,但这并不自动意味着它在筛掉hERG风险这件事上更好,后面还要结合表2和实验验证一起看。 Model 1与Model 2性能对比 靶点 Model 1最佳算法 Model 1 AUC-ROC Model 2最佳算法 Model 2 AUC-ROC ADRB2激动剂 SVM 0.93±0.03 SVM 0.91±0.07 ADRB2拮抗剂 SVM 0.92±0.02 SVM 0.96±0.02 CHRM1激动剂 NB 0.84±0.04 SVM 0.89±0.04 CHRM1拮抗剂 RF 0.94±0.01 SVM 0.96±0.01 DRD2激动剂 SVM 0.88±0.03 SVM 0.90±0.03 DRD2拮抗剂 SVM 0.92±0.02 SVM 0.94±0.03 HTR2A激动剂 SVM 0.84±0.03 SVM 0.86±0.01 HTR2A拮抗剂 SVM 0.92±0.01 SVM 0.94±0.02 hERG抑制剂 SVM 0.91±0.01 NA NA AUC-ROC结果表明大多数模型表现良好,至少有一种ML方法在每个GPCR靶点上达到AUC-ROC>AUC-ROC阈值0.84,在预测hERG抑制剂时达到AUC-ROC=AUC-ROC值0.90 GPCR的AUC-ROC值范围为AUC-ROC下限0.70至AUC-ROC上限0.94,hERG抑制剂的AUC-ROC值范围为AUC-ROC下限0.81至AUC-ROC上限0.91 SVM在大多数GPCR和hERG分类任务中表现最佳,表明其在处理高维分子描述符方面的优势 模型稳定性:10次迭代的性能指标(表S1)显示高度一致性,支持模型达到稳定性能。BAC和MCC的最优值遵循与AUC-ROC相同的趋势,即当AUC-ROC值较大时,BAC和MCC也显示较大值。 骨架拆分验证 为了评估结构泛化能力,研究使用Bemis-Murcko骨架拆分评估了RF和SVM模型。如预期的那样,基于骨架的分区降低了大多数靶点的AUC,反映了预测新型化学类型活性的难度。 CHRM1激动剂和HTR2A拮抗剂观察到最大的下降,可能是由于这些靶点的活性化合物结构多样性有限,限制了骨架特定特征的可转移性。 相比之下,包括ADRB2和CHRM1拮抗剂以及DRD2激动剂/拮抗剂在内的几个靶点的模型保持了相对较高的AUC(AUC下限0.80至AUC上限0.89),表明更一致的结构-活性关系。 总体而言,骨架拆分分析表明,虽然在严格的骨架分离下性能有所下降,但模型对多个GPCR靶点和hERG抑制保留了有意义的预测能力。 结构冗余评估:在使用Tanimoto系数评估LOPAC外部验证集与训练数据之间的结构冗余后,发现630个LOPAC化合物的Tanimoto系数为1,表明可能是重复化合物。这些高相似性化合物可能会高估外部验证性能,因此研究在计算PPV时排除了这些化合物。 外部验证结果 使用LOPAC数据集(Library of Pharmacologically Active Compounds)作为外部验证集评估了在Tox21 10K数据上训练的模型性能。 表1:基于LOPAC实验的两种建模流程外部验证结果 GPCR Model 1最佳算法 Model 1 PPV范围 Model 2最佳算法 Model 2 PPV范围 CHRM1激动剂 SVM 0.41-1.00 SVM 0.47-1.00 CHRM1拮抗剂 SVM 0.65-0.95 SVM 0.64-0.94 HTR2A激动剂 XGB 0.65-0.90 XGB 0.60-0.90 DRD2拮抗剂 SVM 0.74-0.90 SVM 0.73-0.86 ADRB2拮抗剂 RF 0.58-0.81 RF 0.53-0.76 DRD2激动剂 XGB 0.32-0.69 SVM 0.30-0.73 ADRB2激动剂 SVM 0.54-0.64 RF 0.51-0.68 HTR2A拮抗剂 RF 0.14-0.20 RF 0.14-0.23 hERG抑制剂 SVM 0.93 NA NA 外部验证显示大多数模型表现良好,至少有一种ML方法在每个GPCR靶点上达到PPV>PPV阈值0.64(Model 1)或PPV>PPV阈值0.68(Model 2)。SVM在识别hERG抑制剂方面表现突出,Model 1的SVM达到PPV为0.93。值得注意的是,由于原始LOPAC集合中只有5个HTR2A拮抗剂,研究添加了49个经验证的其他活性物质使总数达到54个,产生了更可靠的PPV。 表2:GPCR激动剂与拮抗剂的平均hERG抑制效力 原文的表2比较了不同靶点、不同模式下化合物的平均hERG抑制强度,以 -LogAC50 表示。这个表很关键,因为它回答的不是谁的分类分数更高,而是模型挑出来的分子到底是不是更不容易打到hERG。 靶点 模式 Active Inactive Model 1 active Model 1 active(hERG-inactive only) Model 2 active ADRB2 激动剂 4.32 ± 0.54 4.14 ± 1.00 4.17 ± 0.35 4.12 ± 0.31 3.61 ± 1.99 ADRB2 拮抗剂 4.63 ± 0.63 4.07 ± 1.00 4.73 ± 0.80 4.16 ± 0.42 4.75 ± 0.88 CHRM1 激动剂 4.24 ± 1.09 4.15 ± 0.96 4.27 ± 0.54 4.00 ± 0.00 4.24 ± 0.51 CHRM1 拮抗剂 4.58 ± 0.82 4.03 ± 0.98 4.79 ± 0.66 4.08 ± 0.27 4.65 ± 0.68 DRD2 激动剂 4.31 ± 0.41 4.15 ± 1.00 4.35 ± 0.40 4.17 ± 0.30 4.33 ± 0.40 DRD2 拮抗剂 4.36 ± 1.37 4.05 ± 0.65 4.93 ± 0.75 4.20 ± 0.41 4.92 ± 0.81 HTR2A 激动剂 4.44 ± 1.05 4.06 ± 0.92 4.39 ± 0.51 4.15 ± 0.29 4.53 ± 0.61 HTR2A 拮抗剂 4.32 ± 0.73 4.16 ± 0.97 4.68 ± 0.89 4.17 ± 0.92 4.20 ± 0.75 这张表支持了文中的一个重要判断:GPCR活性化合物,尤其是拮抗剂,平均来看往往伴随更强的hERG抑制;而在Model 1中先用hERG模型做排除,通常能把预测命中的hERG抑制强度再往下压一截。换句话说,Model 2在若干分类指标上略占优,但Model 1在先识别活性、再剔除hERG风险这条路线下,对降低hERG负担更直接。 实验验证结果 图3:模型预测的CHRM1激动剂/拮抗剂的实验验证 图3A-C:代表性强效CHRM1激动剂的结构和浓度-响应曲线,绿色曲线表示CHRM1活性,红色曲线表示hERG活性 图3D-F:代表性CHRM1拮抗剂的结构和浓度-响应曲线,绿色曲线表示CHRM1活性,红色曲线表示hERG活性 图3是全文最重要的落地证据。前三个例子显示,模型不仅能找到CHRM1激动剂,而且这些化合物的绿色曲线与红色曲线明显分开,说明CHRM1活性先出现而hERG作用较弱。后三个拮抗剂例子也传达同样的信息,即真正值得继续推进的,不只是有活性,而是活性与hERG风险之间有窗口。 CHRM1激动剂验证 指标 结果 第一轮测试数量 382个预测CHRM1激动剂 确认为活性 274个 PPV 阳性预测值0.72 强效激动剂 103个,$\mathrm{EC50} < 10~\mu\mathrm{M}$ 代表化合物1 NCGC00642171-01,$\mathrm{EC50} = 1.06 \pm 0.10~\mu\mathrm{M}$ 代表化合物2 NCGC00525960-01,$\mathrm{EC50} = 1.68 \pm 0.50~\mu\mathrm{M}$ 代表化合物3 NCGC00657555-01,$\mathrm{EC50} = 4.21 \pm 1.23~\mu\mathrm{M}$ 这部分结果说明,模型在激动剂方向上的主要价值是把极低的原始命中率显著拉高,并且挑出了一批后续值得进入确认实验的候选。 CHRM1拮抗剂验证 指标 结果 第一轮测试数量 1037个预测CHRM1拮抗剂 确认为活性 945个 PPV 阳性预测值0.91 确认后仍活跃且无显著hERG抑制 66个 强效抑制 34个化合物,$\mathrm{IC50} < 5~\mu\mathrm{M}$ 更强一档 10个,$\mathrm{IC50} < 1~\mu\mathrm{M}$ 已知CHRM1拮抗剂 6个 hERG例外 riboflavin tetrabutyrate 与 NCGC00449480 拮抗剂结果比激动剂更亮眼,尤其体现在PPV上。这也和前面的数据分布一致,即CHRM1拮抗剂数据集本身更大、更容易学到稳定的结构信号。 hERG选择性预测性能 使用阴性预测值(NPV)评估时,TN指在hERG实验中未显示抑制,或hERG抑制效力至少比CHRM1活性低3倍的化合物;FN指以与CHRM1活性相似或更高效力抑制hERG的化合物。总体而言,模型预测化合物在hERG抑制实验中的命中率为命中率18.4%,对应hERG模型的NPV为阴性预测值81.6%。这个结果不能理解成“几乎没有hERG风险”,但足以说明它能把原始化合物库中大量潜在hERG抑制剂预先筛掉。 关键结论与批判性总结 主要贡献 本研究把qHTS数据、QSAR建模、外部验证和后续实验确认串成了一条完整流程。通过比较Model 1与Model 2,作者表明活性预测和hERG风险控制可以被同时纳入同一个筛选框架。对约36万个化合物的虚拟筛选及CHRM1实验证明,这套流程确实能提高命中率,并在一定程度上降低hERG相关风险。实验验证结果显示,ML模型可用于识别具有最小hERG抑制的潜在GPCR药物,模型在识别具有最小hERG抑制的新GPCR靶向化合物方面表现良好。这些模型预测的GPCR靶向化合物为实验测试和进一步开发为药物先导化合物提供了优先级排序的候选列表,为开发更安全的GPCR靶向疗法提供了框架,强调了平衡疗效和心脏安全性的策略需求。 方法学优势 双模型策略:Model 1提供了GPCR与hERG的独立预测,Model 2则把去除hERG抑制剂这一步提前到了建模阶段,两者侧重点不同。根据模型去除hERG抑制剂能力的评估,分别为GPCR靶点和hERG构建的独立模型在去除hERG抑制剂方面比从训练数据中预先排除hERG抑制剂的模型更有效 共识模型:四种ML算法联合决策,减少了单一模型偶然命中的影响。与CardioGenAI和CToxPred2等先进hERG责任框架相比,本研究的分类模型(特别是XGB和SVM)表现出更高的特异性(特异性范围0.98-0.99)和更强的平衡准确率(XGB=平衡准确率0.77,SVM=平衡准确率0.75) 实验闭环:不是停留在交叉验证或外部验证,而是进一步做了CHRM1与hERG实验确认,发现了多个具有微摩尔级活性的新型CHRM1激动剂和拮抗剂,且大多数CHRM1激动剂和拮抗剂对hERG抑制的影响较小(hERG实验中IC50>IC50阈值6.2μM) 可解释比较:不仅比较分类指标,还用表2直接比较了命中化合物的hERG抑制强度,为模型选择提供了定量依据 局限性 仅验证CHRM1:由于资源限制,研究仅对CHRM1预测化合物进行实验验证,其他GPCR模型(ADRB2、DRD2和HTR2A)的实验验证性能可能不同,且一些预测为非活性的化合物可能实际上是活性的(即假阴性) 体外实验依赖性:研究仅应用了一种体外实验方法来生成GPCR靶点和hERG的数据以训练和测试模型,这些实验本身存在假阳性和假阴性率,模型质量因此依赖于实验的技术和生物学可靠性。例如,CHRM1激动剂模式实验的确认率相对较低 单一心脏毒性终点:研究仅考虑了hERG依赖性心脏毒性,未考虑来自其他潜在途径的心脏毒性效应 骨架泛化能力:骨架拆分验证表明模型在预测新型化学类型时性能下降,在某些GPCR靶点(如CHRM1激动剂和HTR2A拮抗剂)观察到最大下降,可能是由于这些靶点的活性化合物结构多样性有限 未来方向 扩展验证范围:对其他GPCR靶点(ADRB2、DRD2、HTR2A)的预测化合物进行实验验证,评估模型在不同靶点上的泛化能力 多目标优化:探索同时考虑GPCR活性、hERG抑制与其他ADMET性质的多目标筛选策略,优化hERG排除阈值以适应不同GPCR靶点和项目阶段的风险容忍度 数据来源多样化:尝试更丰富的分子表示方法和更广的训练数据来源,提升模型对新骨架的外推能力 多心脏毒性终点整合:除了hERG依赖性心脏毒性外,还应考虑来自其他潜在途径的心脏毒性效应,构建更全面的心脏安全性预测框架
Machine Learning & AI
· 2026-03-22
Deep Learning破解双功能抗菌肽设计:DeepQSAR模型的应用与突破
Deep Learning破解双功能抗菌肽设计:DeepQSAR模型的应用与突破 本文信息 标题:The Use of DeepQSAR Models for the Discovery of Peptides with Enhanced Antimicrobial and Antibiofilm Potential 作者:Jiaying You, Hazem Mslati, Evan F. Haney, Noushin Akhoundsadegh, Robert E.W. Hancock, Artem Cherkasov 发表时间:2025年 单位:加拿大不列颠哥伦比亚大学(UBC)、渥太华大学,加拿大 引用格式:You, J., Mslati, H., Haney, E. F., Akhoundsadegh, N., Hancock, R. E. W., & Cherkasov, A. (2025). The use of DeepQSAR models for the discovery of peptides with enhanced antimicrobial and antibiofilm potential. Journal of Chemical Information and Modeling, https://doi.org/10.1021/acs.jcim.5c02138 源代码:https://github.com/chill-bear/peptides(包含数据预处理脚本、模型训练代码和图表生成脚本) 摘要 针对抗生素耐药性的全球危机,抗菌肽(AMPs)因其独特的直接杀菌机制和较低的抗性风险而被视为下一代治疗药物。然而,同时预测抗菌和抗生物膜活性的计算方法仍然匮乏。本研究开发了一种新型Deep QSAR框架,将前馈神经网络(用于定量预测生物膜抑制效率)与递归神经网络(用于二分类抗菌活性预测)相结合,通过迁移学习实现高效的多目标肽设计。模型在5折交叉验证中达到90%的准确率,准确度和召回率分别为0.90和0.88。基于模型预测和聚类分析,我们合成并验证了100个设计肽,其中44个显示出优于参照肽IDR-1018的抗生物膜活性,31个表现出更强的抗菌活性,29个实现了两种活性的同步增强。分子动力学(MD)模拟揭示了这些肽通过强而选择性的细菌膜结合机制(特别是多价的赖氨酸/精氨酸-磷脂酸酯相互作用)来实现其效能的。 核心结论 首次实现同时预测:开发了第一个能够同步预测抗菌和抗生物膜活性的Deep QSAR模型,突破了传统单一功能预测工具的局限 显著的性能优势:90%准确率(5折CV)、精确度0.90、召回率0.88,超过现有Macrel、AI4AMP和DBAASP等预测工具。 实验验证的成功率高:100个设计肽中29个实现双功能增强,展现出模型的实用价值。 最强肽的10倍增强:顶级肽MVLRIKLRLKIR对生物膜的IC50仅为0.147 μM,较参照肽(1.417 μM)低近10倍。 机制清晰:MD模拟证实膜结合和选择性是关键驱动因素,为进一步优化提供了理论基础。 背景 抗生素耐药性(AMR)已成为全球公共卫生危机。过度使用和滥用抗生素导致环境污染加剧,迫使微生物产生防御机制。这不仅削弱了现代医学的治疗效果,还增加了医疗成本和感染死亡率。然而,小分子抗生素的传统开发模式面临瓶颈:新药研发周期长、成本高、成功率低,而且耐药菌株快速进化。 抗菌肽(AMPs)是一类天然防御性蛋白质,广泛分布于细菌、植物、真菌和动物中。与传统小分子抗生素不同,AMPs通过直接破坏细菌膜(如pore formation、carpet机制)和诱导细胞内应激反应来杀灭病原体,这种机制导致耐药菌株发展的可能性大大降低。此外,AMPs还展现出对生物膜的抑制活性,这对治疗慢性感染和医疗器械相关感染至关重要。 然而,传统AMP发现仍依赖于高成本的高通量实验筛选和试错法。虽然机器学习(ML)和深度学习(DL)技术在近年来加速了肽设计过程,但现有的计算工具多专注于单一功能预测(通常是抗菌活性),而对生物膜抑制的预测能力有限。这导致发现同时具有两种功能的候选肽变得极其困难,阻碍了下一代治疗药物的开发。 关键科学问题 本研究针对以下核心问题进行了深入探索: 问题一:如何在同一模型框架中预测多个端点的AMP活性? 传统QSAR模型通常采用单一预测目标(如抗菌活性),基于简单的分子描述符或SMILES编码。而肽的序列和功能的多样性使得多目标预测成为独特挑战——需要模型既能捕捉序列模式信息,又能准确回归生物膜抑制的定量数据。 问题二:如何有效利用异质数据源进行转移学习? 本研究整合了自建的抗生物膜活性数据库(约700个肽,3000个数据点)和三个大规模公开AMP数据库(DRAMP、AI4AMP、DBAASP,共52000余条目)。这些数据来源差异大、标注方式不一、样本分布不均,如何在保证泛化性的同时充分利用这些信息是关键。 问题三:设计的肽能否真正优于参照标准? IDR-1018作为well-characterized的宿主防御肽,已被证明具有广谱抗生物膜活性。新设计肽需要通过实验验证来证明其优越性,这要求模型不仅预测准确,还需识别那些未被充分探索但具有高潜力的序列空间区域。 创新点 首个多功能Deep QSAR框架:融合前馈网络(定量)和RNN(分类),通过迁移学习实现抗菌和抗生物膜活性的同步预测,打破了传统单功能预测工具的局限 创新的双模型整合策略:Model 1的数值输出(生物膜IC50预测)直接馈入Model 2作为额外特征,增强了RNN对序列的理解,避免了简单的模型stacking 大规模数据融合:自建in-house数据库与DRAMP、AI4AMP、DBAASP三大公开库的整合,构建了迄今最全面的AMP训练集,提升了泛化能力 实验验证与机制解析的结合:不仅进行体外活性测试(抗菌、抗生物膜、溶血、细胞毒性),还通过微秒级MD模拟精准解析设计肽的膜相互作用,建立了序列-结构-活性的完整链条 研究内容 核心方法:DeepQSAR双模型框架 为了实现同时预测抗菌和抗生物膜活性,该研究设计了一个创新的两阶段深度学习框架。 graph TB A["肽序列<br/>MVLRIKLRLKIR"] --> B["One-hot编码<br/>每个氨基酸→20维向量"] subgraph "Model 1: 定量预测" B --> C["前馈神经网络<br/>FFNN"] C --> D["In-house数据训练<br/>约700肽, IC50数值"] D --> E["生物膜IC50输出<br/>定量预测"] end subgraph "Model 2: 分类预测" B --> F["递归神经网络<br/>BiLSTM"] F --> G["迁移学习<br/>Model 1权重初始化"] G --> H["公开数据集训练<br/>DRAMP、AI4AMP、DBAASP<br/>52000余个肽, 标签"] H --> I["抗菌活性分类<br/>二分类输出"] end E --> J["特征融合<br/>IC50预测 + 序列模式"] I --> J J --> K["最终预测<br/>双功能评分"] K --> L["聚类与筛选<br/>选择top肽合成验证"] 方法详述: 1.数据准备与编码: 使用One-hot编码将20种标准氨基酸转换为长度为20的二进制向量,保留了序列信息的顺序性。 自建in-house数据集由Kinexus生物信息公司合成的肽阵列组成,约700个唯一肽,每个肽测定了抗MRSA生物膜的IC50值(共3000个数据点)。 结合DRAMP(22259肽)、AI4AMP(10716正例+10718负例)和DBAASP(19751活性肽),构建了超过52000条目的训练集。 2.Model 1—前馈神经网络(定量预测): 基于in-house数据集构建,目标是学习肽序列和生物膜抑制IC50的定量关系 输入为One-hot编码的肽序列,通过多层前馈网络处理,直接输出IC50数值预测 这一模块为后续的RNN提供了生物膜抑制的数值信息基础 3.Model 2—递归神经网络(分类预测): 采用双向LSTM(BiLSTM)架构,针对抗菌活性进行二分类(活性/非活性) 关键创新是迁移学习:将Model 1的训练权重初始化到Model 2 使RNN能够继承关于肽序列和生物膜相互作用的知识,学习序列特定的抗菌模式 融合了定量的生物膜抑制信息,实现多维度特征学习 4.整合与特征融合: 将Model 1的IC50预测输出与Model 2的RNN架构级联(concatenate) 使分类器能够利用数值洞察和序列模式来做出更准确的预测。 设计优势:在架构层面实现了信息流的有机整合,比简单的模型融合更有优势。 数据集与实验方法 使用了四个主要数据源: In-house数据:约700个肽,3000个抗MRSA生物膜IC50测定数据 DRAMP:22259个肽,包含综合的已知AMP序列 AI4AMP:平衡数据集,10716个正例加10718个负例 DBAASP:19751个已实验验证的活性肽 肽的合成在芹菜素膜阵列上进行(Kinexus),通过如下步骤测定活性: 甲氧西林耐药金葡萄球菌(MRSA)用作检验菌株 肽浓度范围1-256 μg/mL,测定OD600(生长)和水晶紫吸收(生物膜) 使用非线性回归拟合IC50值(50%抑制浓度) 模型性能评估 图3:Model 1前馈神经网络的训练过程 前馈网络的平均绝对误差(MAE)和损失函数都随迭代次数逐步下降,最终在验证集上MAE约1.5,表明模型能够准确预测生物膜IC50的量级。训练和验证曲线显示稳定收敛,未出现过拟合现象。 图4:Model 2递归神经网络的分类性能 精确度-召回曲线(左)和ROC曲线(右)显示模型在不同阈值下都保持90%以上的精确度,同时维持88%的召回率。ROC曲线的AUC接近1.0,说明模型具有优异的区分活性和非活性肽的能力。 定量评估结果为: 精确度(正样本):0.90 召回率(正样本):0.88 F1得分:0.89(两个类都>0.88,说明性能均衡) 这些指标远优于现有工具(见附录对Macrel、AI4AMP和DBAASP的对比)。 高通量筛选与设计肽的验证 筛选流程: 从UniProt数据库中提取了20417个已审核的人类蛋白序列 过滤掉长度<100氨基酸的蛋白 使用滑动窗口方法(每次移动1个位置)系统性地分割成12-mer肽片段 对约50000个候选肽进行了预测,筛选出预测评分最高的100个 按照序列相似性进行层级聚类,从每个簇中选择最高评分肽用于化学合成和生物检验 图5:设计肽与训练肽的IC50对比 使用小提琴图展示了设计肽和训练肽在抗生物膜和抗菌两个维度的IC50分布。中位IC50值为: 活性类型 训练肽(μM) 设计肽(μM) 抗生物膜 1.59 0.91 抗菌(浮游) 1.46 1.42 设计肽的生物膜IC50显著低于训练肽,表明模型成功识别并优化了生物膜抑制特性。 图6:合成肽的有效性分类 100个设计肽按照相对于参照肽IDR-1018的表现分类: 44肽:抗生物膜活性更强 2肽:仅抗菌更强 25肽:两者都改善但幅度不同 29肽:两种活性都明显优于对照——这是最有价值的候选 顶级肽的表征 表1:Top 5双功能肽(抗生物膜与抗菌均优) 肽ID 序列 抗生物膜IC50(μM) 抗菌IC50(μM) 10 WKKKGRMRWKWI 0.27 0.74 20 LKIKVHIYRMKR 0.35 1.07 99 MLIRVRKLWRIL 0.24 0.70 40 RARGRKRLVVTI 0.30 1.18 86 RALKKIIKRLCR 0.38 0.70 IDR-1018(对照) VRLIVAVRIWRR 1.42 1.73 最强肽(ID 105, MVLRIKLRLKIR)在抗生物膜上达到0.147 μM,约为IDR-1018的1/10,这代表了迄今最强的AMP生物膜抑制活性之一。其抗菌IC50为1.29 μM,也优于对照的1.73 μM。 图7:阵列肽生物膜vs MRSA活性 该图展示了Top 5肽及对照肽在肽阵列上的生物膜和浮游菌抗性活性曲线。六个面板分别对应肽ID 10、20、99、40、86和105(对照为IDR-1018),每个肽的剂量-反应曲线清晰显示了其多维度效能,进一步验证了设计肽相比对照的改进。 安全性评估 为评估毒性风险,对三个代表肽(J20、J28、J39)进行了溶血和PBMC细胞毒性测定。结果表明: 溶血IC50:全部>250 μg/mL,显示对红细胞的膜破坏极小 PBMC细胞毒性:J28、J39的IC50 >250 μg/mL;J20为166 μg/mL 治疗窗口:生物膜IC50(MBIC)为1-4 μg/mL,远低于毒性阈值,提供了60-250倍的安全边际 这表明设计肽具有良好的生物相容性,适合进一步的临床前开发。 分子动力学揭示作用机制 通过微秒级MD模拟(GROMACS + MARTINI 3粗粒化力场),对43个设计肽在三种膜系统(革兰氏阳性菌模型、革兰氏阴性菌模型、哺乳动物细胞对照)中的相互作用进行了表征。 图8:MD模拟结果——肽-膜相互作用、驻留、选择性和构效关系 A子图 - 时间分辨赖氨酸/精氨酸-膜接触: 所有肽在50-100 ns内建立与膜的多价接触,然后维持高位 抗浮游设计肽:接触数最高(平均15.4,峰值16.7) 双功能肽:中间水平(约13.0) 抗生物膜肽:较低但稳定(约10.8) IDR-1018对照:接近抗生物膜肽(11-12) 非活性肽:无接触(缺乏赖氨酸/精氨酸) B子图 - 磷酸头基团接触密度分布: 磷酸头基团接触密度定义为肽与膜磷脂头基团(PO4)在0.5 nm范围内的接触数,反映肽与膜表面的结合密集程度: 抗浮游菌肽和双功能肽:峰值约3.3 抗生物膜肽:峰值约2.3 IDR-1018:约2.3(与抗生物膜类相同) 非活性肽:仅0.8(极少接触) 设计肽与膜表面的多价磷酸结合密度远高于对照肽,表明肽通过多个精氨酸/赖氨酸残基同时结合多个磷酸基团,形成稳定的多价网络结构,这是膜破坏和细胞溶解的前提条件。 C子图 - 肽-膜中面距离热力图: 热力图显示肽在1微秒模拟过程中与膜的轴向距离演变。根据原文,使用GP膜(革兰氏阳性,用于评估浮游菌杀伤)和GN膜(革兰氏阴性,用于评估生物膜抑制): 抗浮游菌肽和双功能肽:在GP膜上保持浅层驻留(z值约-0.5至0 nm),全程稳定 抗生物膜肽:在GN膜上保持近表层驻留(z值约-0.5至0 nm),持久不变 IDR-1018:界面驻留但波动更大,不如设计肽稳定 非活性肽:远离膜(z值小于-3 nm),无实质接触 D子图 - 选择性评估(细菌膜 vs 哺乳动物膜): 设计肽(所有类):接触数差(Δ)均值约30 contacts(相对于哺乳动物细胞膜),分布集中 IDR-1018:类似正偏移(25-30范围) 非活性肽:接近零(无选择性) E子图 - 构效关系(插入深度vs活性): 肽膜插入深度与活性的相关性分化明显: 抗菌活性(浮游,革兰氏阳性): Spearman相关:ρ = 0.69, p = 0.0045(显著正相关) 趋势:浅层插入与低IC50(高活性)强烈关联 解释:保持在浅表的肽能更有效地破坏膜结构,形成孔隙或地毯溶解;深度插入反而降低活性 抗生物膜活性(革兰氏阴性): 相关性:无显著相关(p > 0.05) 含义:生物膜抑制机制不依赖于膜插入深度,可能依赖于膜表面捕获后的胞内信号干扰(如ppGpp、quorum sensing) Q&A Q1: 为什么One-hot编码而不用其他肽特征(如BLOSUM矩阵、物化性质)? A1: One-hot编码保留了序列的精确顺序信息和完整的氨基酸恒等性,这对RNN学习局部和全局序列模式至关重要。物化性质或BLOSUM会损失肽的某些特异性特征(如某个Cys位置的disulfide潜力)。此外,One-hot编码与循环网络的设计在概念上更贴切——RNN本身就是为处理离散序列而优化的。 Q2: Model 1和Model 2之间的迁移学习具体如何工作? A2: Model 1在in-house抗生物膜数据集上训练,学习了肽序列到IC50(数值)的映射。其中间层权重编码了肽的生物膜亲和力。Model 2初始化时直接复制这些权重到BiLSTM的嵌入层,使RNN一开始就知道哪些序列特征与膜相互作用相关。后续在大型AMP分类数据集上微调时,RNN保留了这些初始化的特征,同时学习抗菌活性的额外模式。这比随机初始化快速得多,也减少了过拟合的风险。 Q3: 为什么选择12-mer作为设计肽的长度? A3: 12氨基酸是最小可行的功能肽长度(short peptides),足以形成α-螺旋或其他二级结构,但避免了合成和成本的复杂性。UniProt滑动窗口方法系统性地生成了大量候选,而12-mer的长度也是文献中well-characterized肽(如IDR系列)的标准。这样既保证了生物学意义,也便于后续的优化。 Q4: 设计肽对其他常见致病菌(如绿脓杆菌、肠杆菌)的活性如何? A4: 论文中仅报告了对MRSA的测定数据(革兰氏阳性)。对广谱活性的验证(包括革兰氏阴性菌)计划在后续研究中进行。MD模拟显示肽在革兰氏阴性模型膜上也有强劲的结合,但体外验证仍是必要的——这也是论文Discussion中强调的局限性。 关键结论与批判性总结 研究意义与影响 开创性的多目标预测框架:首次实现在单一模型中同时预测抗菌和抗生物膜活性,为多功能AMP设计树立了新范式 高实用性的设计管道:从50000个候选肽到100个合成肽,再到29个双功能增强肽,展现了29%的实现率,远超随机合成 强有力的实验验证:不仅测定了生物活性,还进行了毒性评估和分子动力学模拟,建立了序列-结构-活性的完整理解 开源资源分享:代码、数据和模型已上传GitHub,便于学术界复现和扩展 存在的局限性 单一菌株验证:实验仅在MRSA上进行,对其他常见致病菌(绿脓杆菌、鲍曼不动杆菌等)的广谱活性需进一步验证 体内模型缺失:所有活性数据来自体外测定(肽阵列),动物模型和临床相关性评估尚待进行 机制理解仍需深化:虽然MD模拟提供了膜相互作用的线索,但关于肽的具体杀菌模式(是否形成孔隙、地毯机制还是其他)仍需要补充生物物理学实验 长期稳定性未评估:肽的血清稳定性、给药形式和体内代谢还没有系统研究 未来研究方向 扩展菌种覆盖:针对多重耐药菌(MDR)、泛耐药菌(XDR)进行活性测定,包括临床分离株 动物模型验证:利用小鼠感染模型评估体内疗效和毒性,为临床前开发奠定基础 结构优化循环:基于MD洞察,进行理性的点突变和截断,进一步提升特异性和效能 AI模型迭代:整合更多数据源(如微生物组数据、宿主防御肽文献),开发下一代多参数预测模型
Machine Learning & AI
· 2025-11-09
DeepQSAR抗菌肽发现——技术细节与扩展数据
DeepQSAR抗菌肽发现——技术细节与扩展数据 完整数据集描述 In-house抗生物膜数据库 约700个唯一肽(多数为12-16氨基酸),由Kinexus生物信息公司通过肽阵列合成(SPOT-array technology)。每个肽针对MRSA进行了2折串联稀释测定,产生了3000个IC50数据点。 数据特征: IC50范围:0.09-50 μM(中位数~1.5 μM) 肽长度分布:8-18氨基酸为主,12-14mer最多 化学修饰:C端酰化(标准AMP格式),某些肽含有非标准氨基酸如Nle(仲亮氨酸)、Trp衍生物 DRAMP 3.0 (Database of Antimicrobial Peptides) 包含:22259肽条目 来源:已发表文献中已知的AMP,涵盖细菌、真菌、植物、昆虫、哺乳动物来源 标注:二进制(活性/非活性),基于文献报道的MIC或IC50阈值 优势:高覆盖度,包括多种菌种的活性信息(不仅限MRSA) 局限:某些条目可能基于定性描述而非精确数值 AI4AMP (Antimicrobial Peptide Predictor) 包含:平衡数据集,10716正例(已知活性AMP) + 10718负例(非AMP序列) 来源:公开AMP数据库与生成的非AMP背景 特点:经过特征工程优化(physicochemical property encoding) 用途:在本研究中主要用于验证和外部基准测试 性能(来自原始论文):精确度~90%,泛化性好 DBAASP v3 (Database of Antimicrobial Activity and Structure of Peptides) 包含:19751活性肽,附带实验验证的结构和活性数据 数据质量:高,仅收录已发表、经实验验证的肽 附加信息:包含部分肽的3D结构、膜交互描述符(如hydrophobic moment、charge distribution) 覆盖范围:广谱菌种(需要标准化处理) 数据集组合与预处理 四个数据源合并后,采用如下预处理步骤: 去重:基于精确序列匹配移除重复肽 长度过滤:保留8-20氨基酸,去除超短(<8aa)或超长(>20aa)肽,使分布更均匀 编码规范化:将所有非标准氨基酸(如Nle、Orn)映射到最相近的标准氨基酸(Leu、Lys) 标签一致化:对于在多个库中重复出现的肽,采用多数票法决定标签;如信息矛盾则排除 数据平衡:对于分类任务(Model 2),使用SMOTE或加权损失函数处理类不平衡 最终数据集规模:约52000个条目(去重后),其中正例(活性AMP)约占55% 详细方法学 Peptide Clustering算法 为减少合成肽的冗余性并保证序列空间的多样性覆盖,使用了层级聚类(Hierarchical Clustering): 相似性计算:对所有候选肽对进行全局序列比对(Needleman-Wunsch算法),计算相似度矩阵 聚类方法:AgglomerativeClustering(sklearn),使用欧式距离和完全链接(complete linkage) 聚类数:设置为100,对应最终的合成肽数量 代表选择:从每个簇中选择模型预测评分(combined score)最高的肽 优势:确保了100个合成肽均匀分布在5万个候选肽的序列空间中,最大化了发现新功能肽的概率 分子动力学模拟参数 软件和力场: MD引擎:GROMACS 2021.5 粗粒化力场:MARTINI 3.0(适合微秒级长模拟) 初始结构制备:α-螺旋(PeptideBuilder)→ martinize2转换 膜系统构建: 革兰氏阳性菌(GP)膜:POPG:Cardiolipin = 3:1(代表革兰氏阳性菌的外膜) 革兰氏阴性菌(GN)膜:POPE:POPG:Cardiolipin = 6:2:1(代表革兰氏阴性菌的内膜) 哺乳动物对照(MAM):100% POPC(代表人类红细胞膜,用于评估选择性) 模拟条件: 系统尺寸:~15×15×35 nm³ 离子浓度:0.15 M NaCl 温度:323 K(50°C,适合MARTINI) 压力:1 bar(semi-isotropic) 时间步长:20 fs(粗粒化允许) 运行时间:1 μs/复制本,3个复制本/肽/膜(共9 μs/肽) 模拟后分析: 肽-膜接触数(0.5 nm cutoff) Lys/Arg-磷酸基团相互作用(多价结合) 肽中心质量(COM)与膜中面的距离(评估插入深度) RMSD/RMSF(结构稳定性) Spearman相关分析:深度 vs log(IC50),评估插入-活性关系 结果验证:使用MDAnalysis (Python)进行轨迹解析,所有时间序列数据经3个复制本平均后,仅用未平滑数据进行统计(只有图中的类别均值经高斯平滑σ=5) Top 10肽完整列表 Table 1: 最强10个抗生物膜肽 肽ID 序列 抗生物膜IC50(μM) 说明 105 MVLRIKLRLKIR 0.147 最强,约IDR-1018的1/10 39 RGFVRLKKWFNI 0.23 含Trp,可能增强膜插入 99 MLIRVRKLWRIL 0.24 双功能候选(也在抗菌Top 10) 10 WKKKGRMRWKWI 0.27 高Lys密度,强静电结合 59 FRVCYRGICYRK 0.30 含Cys,可能形成disulfide 40 RARGRKRLVVTI 0.30 双功能候选 28 FRVCYRGICYRR 0.35 精氨酸富集,膜结合强 20 LKIKVHIYRMKR 0.35 双功能候选,含疏水残基 86 RALKKIIKRLCR 0.38 双功能候选,平衡疏水-亲水 IDR-1018(对照) VRLIVAVRIWRR 1.42 参照标准 Table 2: 最强10个抗菌(浮游)肽 肽ID 序列 抗菌IC50(μM) 说明 99 MLIRVRKLWRIL 0.70 最强,双功能 86 RALKKIIKRLCR 0.71 双功能,高效率 10 WKKKGRMRWKWI 0.74 双功能 102 VLRIGWILWRIS 0.84 高疏水性 62 RRRAKGRIRLIV 0.89 Arg富集 100 LLILWRKLWILR 1.02 疏水性主导 2 GRMRWKWIKKRI 1.03 基础设计 20 LKIKVHIYRMKR 1.07 双功能 33 GLKSFARVLKKI 1.15 序列多样性 40 RARGRKRLVVTI 1.18 双功能 IDR-1018(对照) VRLIVAVRIWRR 1.73 参照标准 关键观察: 5个肽同时出现在两个Top 10中(ID 10, 20, 40, 86, 99),这些是最有价值的候选 抗生物膜肽倾向于高Lys/Arg密度和Trp含量(增强膜亲和力) 抗菌肽显示更多的疏水残基组合(增强膜插入和破坏能力) 与其他AMP预测工具的对比分析 三种现有工具的性能 本研究在29个实验验证优于IDR-1018的肽上,对比了三个广泛使用的AMP预测工具: Macrel (AMP Mining in Genomes and Metagenomes) 原理:22个物化描述符(电荷、疏水性矩、二级结构倾向等) + 传统ML分类器 结果: 29个验证肽的预测评分范围集中在0.50-0.60区间 接近默认阈值(0.50),导致低区分度 假阴性率高,精确度~50%,召回率同样低 局限:Macrel设计用于基因组/宏基因组挖掘(未知序列背景),对已知AMP数据库的表现不理想 AI4AMP (Antimicrobial Peptide Predictor) 原理:物化性质编码 + 卷积神经网络(CNN) 性能: 在定性上,对多数29个肽给出了高AMP概率评分 但当以IDR-1018的评分作为分类阈值时,精确度和召回率均~50% 混淆矩阵显示该阈值选择不当,导致过多假阳性或假阴性 优点:模型本身性能不错,但对于高活性肽的定量区分有限 DBAASP Predictor 原理:三个膜交互相关描述符(hydrophobic moment、charge density、membrane-depth potential) 结果: 29个肽中,正负预测几乎均分(接近50:50) 基于这三个特征的区分能力有限 虽然这些描述符在AMP设计中重要,但单独使用不足以预测多功能性 反思:强调了序列-序列相关性(通过RNN捕捉)的重要性,单纯依靠物化特性难以抓住功能差异 DeepQSAR的优势总结 指标 Macrel AI4AMP DBAASP DeepQSAR 精确度 ~50% ~50% ~50% 90% 召回率 低 低-中 低 88% F1得分 <0.5 0.40-0.50 <0.5 0.89 多目标预测 否 否 否 是 泛化性 有限 中等 一般 优异 计算成本 低 中 低 中-高 毒性与安全性数据 溶血活性 三个代表肽(J20、J28、J39)在人红细胞上的溶血测定: 图S1展示的浓度-反应曲线表明: J20 (LKIKVHIYRMKR):IC50 >250 μg/mL(上限未达),极低溶血风险 J28 (FRVCYRGICYRR):IC50 >250 μg/mL J39 (RGFVRLKKWFNI):IC50 >250 μg/mL 解释:即使在256 μg/mL(最高测试浓度),红细胞溶解也<10%,说明对宿主细胞膜的破坏最小。相比之下,许多阳性对照AMP在10-50 μg/mL即表现出明显溶血。 PBMC细胞毒性 外周血单核细胞(PBMCs)对肽的耐受性评估: 数据来自Table S1: | 肽 | PBMC IC50(μg/mL) | 与MBIC的倍数差 | 评价 | |—-|—————-|————-|——| | J20 | 166.1 | 41-166倍 | 中等毒性 | | J28 | >250 | >62.5-250倍 | 低毒性 | | J39 | >250 | >62.5-250倍 | 低毒性 | 最小生物膜抑制浓度(MBIC):1-4 μg/mL(与IC50测定相同条件) 治疗窗口:毒性IC50 / MBIC = 62-250倍,足以支持临床前开发(理想值通常>10倍) PBMC毒性的分化原因: J28/J39高度耐受,可能与其特定的Cys、Tyr组成(可能稳定膜界面而不破坏)有关 J20的中等毒性可能源于其高Lys密度,在高浓度时对人细胞也有一定膜扰动 补充图表详解 Figure S1: 溶血活性曲线 左图为Hemolysis,右图为PBMC Cytotoxicity,横轴肽浓度(log scale, 0.6-256 μg/mL),纵轴为百分比溶解/毒性。三条曲线代表J20(蓝)、J28(红)、J39(绿)。 关键发现:三肽在1-256范围内溶血均<15%,PBMC毒性中J28/J39始终<10%,J20在128-256 μg/mL才明显上升。 Figure S2: Macrel预测分布 柱状图显示29个验证肽的Macrel评分分布。评分集中在0.50-0.60,大多聚集在单一柱子(34.5%),显示低区分度。 Figure S3: AI4AMP概率评分 曲线图显示概率分布,大多肽评分在0.7-1.0(高AMP概率),但相对于IDR-1018基准(虚线)的区分不足。 Figure S4: DBAASP混淆矩阵 左侧混淆矩阵显示,DBAASP的预测与实际结果的吻合度低,正负预测几近等分。 数据获取与复现 所有数据、代码和预训练模型已公开发布在GitHub仓库: 地址: https://github.com/chill-bear/peptides 内容: data/: 原始IC50数据(CSV)、聚类结果、验证肽序列 models/: 预训练的Model 1和Model 2权重(HDF5格式) scripts/: One-hot编码、模型训练、超参数调优、图表生成代码(Python) md_simulations/: MD设置文件(.top, .gro, .mdp)、轨迹分析脚本 复现步骤: Clone仓库并安装依赖(TensorFlow, scikit-learn, MDAnalysis等) 运行预处理脚本整合四个数据源 使用提供的超参数训练Model 1和Model 2 对自有候选肽进行预测和聚类 用GROMACS运行MD模拟,使用MDAnalysis脚本分析 注意:MD模拟计算密集,建议使用GPU集群或HPC资源;单肽1 μs的三复制本约需2-4小时(单CPU)。
Machine Learning & AI
· 2025-11-09
<
>
Touch background to close