Python神经网络、随机森林、PCA、SVM、KNN及回归实现ERα拮抗剂、ADMET数据预测如果你做过药物发现相关的计算工作一定会遇到这个非常现实的场景手里拿到一批化合物的SMILES结构式生物活性数据零零散散接下来老板就一句话——“先跑几个机器学习模型看看能不能做个预测”。然后就没了。当你自己动手的时候才会发现这里面每一步都藏着坑。从分子怎么变成数字到特征怎么降维再到模型调参、结果评估任何一个环节处理不好最后跑出来的指标都会好看得离谱但实际预测效果一塌糊涂。我今天想完整复盘一下当初用Python完成的一个ERα拮抗剂活性预测和ADMET属性预测项目把中间所有关键的决策点逐个拆开来讲。项目覆盖了分子描述符预处理、PCA降维、随机森林、SVM、KNN、神经网络分类器以及回归模型的构建与验证数据用的是公开化合物库和ADMET基准数据集。整理出来的东西不仅是代码怎么写的重点是在每个环节“为什么这么做”以及我踩过哪些坑。1. 从理化性质到活性终点ERα拮抗剂预测这件事的本质1.1 为什么ERα拮抗剂预测是经典的QSAR问题ERα雌激素受体α是乳腺癌治疗中一个非常重要的靶点。临床上一线内分泌治疗药物比如他莫昔芬、氟维司群本质上都是通过调控ERα信号通路发挥作用。而ERα拮抗剂的意思就是这类化合物结合到受体上之后不激活下游转录反而把通路“按死”从而抑制雌激素依赖的肿瘤细胞增殖。从计算化学的角度看一个化合物到底有没有ERα拮抗活性取决于它的三维结构、疏水性质、氢键供受体分布能不能和ERα的配体结合口袋产生有效相互作用。这些东西听起来很复杂但在QSAR定量构效关系框架下可以把它抽象成一个标准的有监督学习问题输入分子的结构表征向量描述符/指纹输出拮抗活性标签有/无或者IC50/pIC50连续值所以这个项目的本质不是什么生物机制研究而是把药物化学里积累的实验数据拿过来用机器学习模型去拟合“结构-活性”之间的映射关系。1.2 ADMET预测在管线中的真实定位光预测活性还不够一个化合物就算体外活性再强如果吸收差、代谢快、毒性大也成不了药。ADMET就是五个词的首字母Absorption吸收、Distribution分布、Metabolism代谢、Excretion排泄、Toxicity毒性。ADMET预测的目标是在化合物进入动物实验之前先用计算模型把成药性风险高的分子筛掉一批从而节省大量实验成本。在这个项目里ADMET数据集包含了一系列连续或二分类的终点比如水溶性LogS、血脑屏障透过率BBB、人血浆蛋白结合率PPB、hERG心脏毒性等。每个终点单独建一个模型而建模的管线其实和ERα活性预测高度类似都是先算特征再跑模型。这也是为什么这个项目把所有方法放在一起讲因为本质就是一套pipeline换不同的标签反复用。1.3 结构化梳理整个项目的技术栈和流程在执行层面项目核心依赖库是scikit-learn、pandas、numpy、matplotlib神经网络部分用了Keras/TensorFlow或者PyTorch。整个流程可以分为六大块数据准备收集SMILES结构式清洗、去重、去无机盐、剔除无法解析的结构分子表征用RDKit计算分子描述符或生成Morgan指纹把化学结构变成数值向量特征工程零方差过滤、相关性分析、PCA降维模型构建分别用随机森林、SVM、KNN、神经网络、回归算法训练评估验证交叉验证、独立测试集验证分类问题看AUC/Accuracy/Precision/Recall回归问题看R²和RMSE结果分析特征重要性排序、PCA载荷分析、模型的适用域判断。这套流程所有代码都可以在本地跑通只要有Python环境和RDKit就能复现。下面每一节我会把每个环节的细节和处理逻辑单独展开说清楚。2. 数据端的一步都不能马虎SMILES清洗、描述符计算与训练集划分2.1 从原始数据到可用训练集这件事比建模更耗时我做这个项目的时候数据清洗花掉的时间比后面所有建模加起来都长。数据源一般来自ChEMBL或者文献补充材料但拿到的原始SMILES惨不忍睹有的是混合盐形式比如CN(C)CCOC(c1ccccc1)c1ccccc1.Cl是盐酸盐有的是特殊标记比如[Na]CC(O)O还有的SMILES虽然能解析但电荷状态、立体化学标记不统一。这一步的行业标准做法是用RDKit做标准化我通常这样做from rdkit import Chem from rdkit.Chem import Descriptors, AllChem, SaltRemover remover SaltRemover.SaltRemover() def standardize_smiles(smiles): try: mol Chem.MolFromSmiles(smiles) if mol is None: return None mol remover.StripMol(mol, dontRemoveEverythingTrue) smiles_std Chem.MolToSmiles(mol) return smiles_std except: return None标准化之后还需要去重。同一个化合物在数据库里可能出现多次活性值还不一样处理规则一般是生物学重复取平均值但如果数值差异特别大比如超过一个数量级我会把这一条标记为可疑数据宁可删掉也不强行合并。这里要提醒一个非常容易被忽略的点先切分数据集再做特征选择或降维。在完整数据集上算PCA或者选特征信息会从训练集“泄漏”到测试集里导致模型评估虚高。后面我会专门用一节来讲这个坑。2.2 描述符和指纹分子怎么变成机器能读的向量这是整个项目的起点也是最影响模型性能的一步。分子表征主要有两大流派描述符流派用RDKit计算分子的理化性质像MolWt分子量、LogP脂水分配系数、HBD氢键供体数、HBA氢键受体数、TPSA极性表面积、RotBonds可旋转键数等等。这些描述符有明确的物理化学意义计算快维度低通常一两百个模型可解释性比较强。缺点就是它们只描述整体性质丢掉了局部结构信息。指纹流派把分子结构转成一种bit向量。最常用的是Morgan指纹ECFP这种指纹把每个原子周围的环境半径编码成圆形子结构半径2的时候叫ECFP4半径3的时候叫ECFP6。指纹维度非常高默认2048位但包含的结构碎片信息丰富特别适合活性分类这类问题。from rdkit.Chem import AllChem def compute_fingerprints(smiles_list, radius2, n_bits2048): fps [] valid_indices [] for i, smi in enumerate(smiles_list): mol Chem.MolFromSmiles(smi) if mol is None: continue fp AllChem.GetMorganFingerprintAsBitVect(mol, radiusradius, nBitsn_bits) arr np.zeros((n_bits,), dtypenp.int8) from rdkit.DataStructs import ConvertToNumpyArray ConvertToNumpyArray(fp, arr) fps.append(arr) valid_indices.append(i) return np.array(fps), valid_indices我的经验是分类问题优先用Morgan指纹回归问题优先用描述符。原因在于活性分类通常是结构碎片决定活性有无指纹能捕捉到关键的药效团信息而ADMET的连续值预测比如LogS整体的理化性质起主导作用描述符就很够用。2.3 正负样本怎么构造ERα拮抗剂数据集的特殊难点如果你是做活性预测就会知道数据不平衡是个常态。ChEMBL里ERα拮抗剂的数据虽然不少但如果你把“活性阈值”设定在pIC50 ≥ 6.5约IC50 ≤ 316 nM那活性化合物可能只占总数据的三分之一左右。直接用原始数据训练模型会偏向预测“无活性”因为这样准确率很高但毫无意义。我当时的处理办法是把同一靶点的多种活性终点统一成二分类标签取IC50/EC50/Ki数据换算成pIC50之后阈值以上算“positive”阈值以下算“negative”。这样做的好处是用连续值统一了不同实验来源的量纲避免了只保留IC50导致数据量骤减的问题。如果数据量足够还可以考虑用生成对抗网络或者SMOTE做正样本过采样。但在这个项目里我更推荐保守路线先通过加权类别权重缓解不平衡问题同时用AUC而不是Accuracy作为主评估指标因为AUC对类别不平衡不敏感。后面评估部分我会详细讲。3. PCA降维特征从几千变成几十模型反而更稳的原因3.1 多一个特征不一定更好维数灾难的直觉理解很多人刚开始学机器学习时会有一个错觉觉得特征越多模型学到的信息越多效果应该越好。但如果你真的试过2048位Morgan指纹直接扔进KNN或者神经网络大概率会遇到两种情况训练耗时暴增但效果没提升甚至测试集表现更差。原因就是维数灾难。高维空间里样本之间的距离会变得“又远又均匀”原本在高维空间中代表相似子结构的分子其欧氏距离并没有明显区分度KNN完全失效。而SVM在高维空间中寻找超平面也需要更多样本支持否则就是过拟合。PCA在这里的作用是找到一个低维子空间使得数据在这个空间上的投影方差最大也就是用更少的维度保留尽可能多的信息。这相当于把2048位的高维指纹“压缩”成几十个主成分从而让后续的模型在更稠密的低维空间中学习。3.2 PCA具体怎么操作以及主成分个数怎么选PCA的操作本身非常简单scikit-learn里只需要几行代码。但关键问题是到底保留多少个主成分。我以前见过有人直接固定选50个这也太粗暴了。标准做法是先画出“累计方差解释率曲线”看前n个主成分能解释多少比例的方差。一般药剂学相关项目经验是累计解释比例达到70%-80%的时候信息就有了保障这时候对应的主成分数大概是几十个。from sklearn.decomposition import PCA pca PCA(n_components0.8, svd_solverfull) X_pca pca.fit_transform(X_train) print(f原始特征维度: {X_train.shape[1]}) print(f降维后维度: {X_train.shape[1] if hasattr(pca.n_components_, __len__) else X_pca.shape[1]})这里有个容易踩坑的细节n_components0.8在PCA里表示保留80%的方差但你算出来的降维后维度是根据训练集得到的。如果用完整数据先算PCA再切分训练集和测试集之间就发生了泄漏。正确做法是先对训练集fit_transform再用训练集的PCA参数去transform测试集X_test_pca pca.transform(X_test)至于PCA是否可以提升分类准确率我的经验是看模型类型。对KNN来说PCA几乎总是有帮助因为降维重构出的低维表示把最重要的结构差异放在前几个主成分上距离计算更有区分度。对随机森林来说PCA未必有明显增益因为树模型对特征维数和相关性的容忍度很高有时原始特征反而保留更多局部细节。对SVM和神经网络来说PCA通常能带来训练速度和稳定性上的增益尤其是样本数量不多的时候。3.3 PCA之外低方差过滤和相关性去冗余在主成分分析之前还有一道前置处理就是删除“废话特征”。2048位指纹里有大量位在所有样本中要么全是0、要么全是1这些位没有任何区分能力直接删除。另外还有高度相关的一对特征比如两个描述符相关系数达到0.98以上留一个就行。我用的是scikit-learn里的VarianceThreshold但要注意它默认会删除方差为0的特征而且方差阈值需要自己调。由于指纹是二进制向量一个位出现的频率接近0或接近1时它的方差都接近0所以VarianceThreshold(threshold0.01)这个经验值够用。PCA一般放在方差过滤和标准化之后跑。4. 五种算法同台竞技RF、SVM、KNN、神经网络、回归的建模差异4.1 随机森林为什么它永远是药物预测的baseline之王随机森林在QSAR和ADMET预测中的地位就跟SQL在数据处理中的地位一样是绕不开的。它由多棵决策树集成每棵树在训练时用Bootstrap采样出不同的子集并在每个节点分裂时随机选取一部分特征来寻找最佳分裂。这两重随机性使得单棵树之间的相关性降低集成后的模型方差更小、泛化能力更强。用RandomForestClassifier调参时真正需要认真调的参数其实就几个n_estimators树的棵数。我一般设300到500就收敛了再大只会增加训练时间精度提升微乎其微max_features每次分裂考虑的特征数。对于高维指纹数据设置成sqrt特征数的平方根效果比较稳min_samples_leaf叶节点最少样本数。这个参数是抑制过拟合的关键推荐设成5到10尤其当训练数据只有几百个样本时。还有一个容易被忽略的技巧随机森林输出的是样本属于各类别的平均概率而不是投票的硬标签。在二分类场景下把predict_proba得到的阳性概率拿去做AUC计算远比直接用predict的结果更有区分度。这也是我在项目里把随机森林结果当软分类器用的原因。4.2 SVM核函数选择和C值调优的实战逻辑SVM的基本思想是在特征空间中寻找一个“间隔最大”的超平面来分割数据。对于药物活性数据绝大多数情况下特征和标签之间不是线性关系所以要用核函数把原始特征映射到更高维的空间再找超平面。项目里我用的核函数是RBF径向基函数$$K(x, x) \exp(-\gamma | x - x |^2)$$这里有两个关键超参数正则化系数C和核系数γ。C越大模型越不容忍分类错误训练集上越准但越容易过拟合C越小模型越“宽容”泛化可能更好但欠拟合风险增加。γ控制单个样本的影响半径γ越大决策边界越复杂越容易过拟合。我当时用的是网格搜索加5折交叉验证来找最优组合搜索范围是C∈{1, 10, 50, 100}γ∈{1e-3, 1e-2, 1e-1, 1}。一个细节是SVM对特征的尺度极其敏感RBF核里算的是欧氏距离如果特征一个量级是几百一个量级是0.1大的特征会完全主导。所以用SVM之前一定要做标准化通常用StandardScaler。from sklearn.preprocessing import StandardScaler from sklearn.svm import SVC from sklearn.model_selection import GridSearchCV scaler StandardScaler().fit(X_train) X_train_scaled scaler.transform(X_train) X_test_scaled scaler.transform(X_test) param_grid {C: [1, 10, 50, 100], gamma: [1e-3, 1e-2, 1e-1, 1]} svm_model SVC(kernelrbf, probabilityTrue, class_weightbalanced) grid GridSearchCV(svm_model, param_grid, cv5, scoringroc_auc) grid.fit(X_train_scaled, y_train)SVM在中小规模数据集上表现得非常强悍尤其当样本量在几百到几千时它的泛化能力常常优于随机森林。但如果样本量到了几万甚至更多SVM的训练时间会让人崩溃这时候随机森林或神经网络反而更有优势。4.3 KNN算法最简单但对预处理最敏感KNN的原理一句话就能说清楚新样本的标签由它在特征空间里最近的K个邻居投票决定。它不需要训练过程所有“学习”都发生在预测时。但这个项目里KNN一开始效果很一般后来找到的问题是出在标准化上。Morgan指纹的每一位几乎是0/1二值但如果混入了分子量、LogP这种连续描述符量纲完全不同直接算欧氏距离时连续特征权重大得离谱。解决办法有两种要么只用二值指纹做KNN要么对全特征做标准化。我这里推荐一个更稳妥的方案KNN前的特征统一做MinMaxScaler或者StandardScaler在PCA降维之后效果更好。K值的选择也很有讲究。K太小模型对噪声敏感一个离群点就能改变结果K太大远处不相关样本也参与投票决策边界过于平滑。我一般用GridSearchCV扫K从1到21的奇数配合距离加权weightsdistance效果好过用默认的均匀投票。4.4 神经网络药物数据量下怎么避免过拟合在药物发现领域神经网络这个概念通常指的不是深度学习大模型而是多层感知机MLP。MLP由输入层、若干隐藏层、输出层组成每层包含若干神经元层与层之间用全连接方式相连每个神经元做加权求和后经过非线性激活函数。这个项目的神经网络分类器结构比较简单from tensorflow.keras.models import Sequential from tensorflow.keras.layers import Dense, Dropout, BatchNormalization model Sequential([ Dense(128, activationrelu, input_shape(X_train_pca.shape[1],)), Dropout(0.3), Dense(64, activationrelu), Dropout(0.2), Dense(1, activationsigmoid) ]) model.compile(optimizeradam, lossbinary_crossentropy, metrics[AUC])结构上不需要非常深两到三个隐藏层就足够拟合大部分QSAR数据。这里真正需要关注的是过拟合的问题化合物数据集往往只有几百到几千个样本而神经网络的参数成千上万很容易把训练集背下来。我用的防过拟合策略有三个Dropout随机失活、EarlyStopping早停、以及BatchNormalization。尤其是EarlyStopping设monitorval_loss, patience30验证集loss连续30个epoch不下降就自动停止。另外训练集和验证集的比例在神经网络里建议做到8020甚至9010因为数据太少时验证集会不够稳定。神经网络的优势在于预测概率的分布比较平滑不像随机森林那样输出是离散化的概率区间。缺点是训练不稳定同样的种子和结构几次跑下来AUC可能波动几个百分点所以在最终报告里我会多次重复训练取均值而不拿单次结果当结论。4.5 回归模型ADMET连续值预测不能只跑线性ERα拮抗剂活性预测是分类问题但ADMET预测里很多终点是连续数值。比如LogS是log scale的水溶性数值PPB是0到100的百分比。这时候需要的是回归模型。回归这块项目里我实际尝试了三种线性回归Ridge、随机森林回归和SVR支持向量回归。先说结论Ridge回归通常只能作为下限基准它的预测能力受限于特征和标签之间的线性相关性在ADMET任务上非常弱。随机森林回归因为能捕捉非线性关系基本上是默认选择。SVR在中小样本上表现和随机森林接近但调参更麻烦。回归模型的核心评估指标是R²和RMSE。R²越接近1说明模型解释了大部分方差但要注意R²可以因为数据量小或者分布偏斜而虚高。RMSE是有物理单位的误差比如LogS的RMSE如果是0.5意味着平均预测误差在0.5个log单位换成溶解度就是大约3倍的误差。对于ADMET预测RMSE比R²更能反映模型的实际可用性。用随机森林回归时一个实用技巧是设置oob_scoreTrue即用袋外样本做验证。相当于在训练过程中顺便得到一个几乎无偏的验证指标不需要额外切分数据对样本量小的项目特别友好。5. 评估指标怎么读活性分类和ADMET回归的判定逻辑5.1 分类模型不能只看准确率AUC、Precision、Recall的配合使用如果你做的是不平衡的活性分类比如只有25%的化合物有活性准确率会给你一个很大的错觉。全部预测为无活性就已经有75%准确率但这样的模型没有任何实用价值。真正需要关注的指标组合是这组AUC-ROC曲线下面积它衡量的是模型对正负样本排序能力的优劣随机猜测是0.5完美是1.0。AUC不关心阈值怎么定所以它在类别不平衡时依然稳健。这个项目的ERα二分类模型AUC在0.85以上算初步可用0.90以上才是真正有参考价值的模型。Precision精确率预测为活性的化合物里真正有活性的比例。在药物筛选场景Precision低意味着你送去做实验的化合物大多数是白做的实验成本浪费严重。Recall召回率所有真正有活性的化合物里被模型找出来的比例。Recall低意味着活性化合物被漏掉了这比误报更可惜因为发现一个先导化合物的机会可能就这样溜走了。F1 ScorePrecision和Recall的调和平均值在两者之间做一个平衡。我的经验是活性分类项目里如果只能跟别人汇报一个数报AUC如果是给生物学家做筛选还需要告诉他们Precision和Recall的取舍关系。比如我把分类阈值调高Precision上来了但Recall下去了这时候筛选命中率高但漏掉一些潜在活性分子阈值调低则反过来。最终用哪组阈值取决于下游实验成本和对漏检的容忍度。5.2 交叉验证与独立测试集别让模型“偷看”答案这个项目里我始终坚持一个流程先把数据划分成训练集80%、验证集10%、独立测试集10%。GridSearch和EarlyStopping只允许在训练集验证集上做独立测试集从头到尾只能碰一次等到最终模型确定了才用来评估。交叉验证的价值在于利用有限的样本反复验证模型稳定性。常用的有5折和10折我在这个项目里用的是StratifiedKFold分层抽样确保每一折里正负样本比例和全局一致。有时候会发现一个很有意思的现象某几折AUC非常高某几折明显偏低。这种波动不能只看均值还要看标准差。如果标准差超过0.05说明模型对训练数据的构成太敏感这时候应该回去检查数据是否有重复结构比如同一骨架的相似分子被同时分到了训练和测试里。5.3 适用域问题为什么测试集指标好看但真实预测不准最后一个可能让你在项目汇报时翻车的问题模型的适用域。QSAR模型不是在所有化学空间上都适用。你的训练数据如果主要覆盖了含氮杂环类化合物那模型对含硼化合物或者含碘化合物的预测就属于外推可信度极低。我当时在ADMET回归模型里加了一个简单的适用域判定对每个待预测样本计算它到训练集样本的平均欧氏距离在PCA降维后的空间里。如果平均距离超过训练集样本之间平均距离的某个阈值一般是1.5倍到2倍就在预测结果里标记为“低置信度”。这个方法不复杂但非常实用可以避免把模型用在明显不适用的情况下。6. 复盘那些代码之外的事数据泄漏、不平衡与集成思路6.1 我遇到过的数据泄漏那段痛苦的排查经历这个项目的第一次建模结果异常完美AUC高达0.98。当我满心欢喜地拿独立测试集验证时表现依然很好。但往深了想就觉得不对劲因为ERα拮抗剂数据本身噪声很大不可能有这么好的可分性。排查之后发现问题出在一个特别隐蔽的地方我在特征计算之前先对整个数据集做了PCA然后再切分训练集和测试集。这意味着测试集的信息准确说是特征空间的方差结构已经参与到了训练集的降维变换中。这不是一个小问题它会系统性地高估模型性能。第二次我把PCA放在切分之后再跑AUC立刻降到0.87。这个差距就是数据泄漏造成的。类似的泄漏还可能出现在对完整数据做标准化、做缺失值填充、做特征选择。所有“学习型”的数据处理步骤都必须只基于训练集然后应用到测试集。这是我花了两天才想明白的教训写在这里希望你们不用再踩一遍。6.2 正负样本不平衡的应对不是所有数据增强都合适前面说过ERα拮抗剂数据集天然不平衡。我当时对三种策略做了对比类别权重、SMOTE过采样、随机欠采样。结果是SMOTE在训练集上表现不错但独立测试集上并没有明显优于单纯使用类别权重的方法反而增加了过拟合风险。随机欠采样会损失大量样本也不推荐。最后我选择的是class_weightbalanced。这个方法简单直接让模型在训练时对少数类的样本错误付出更大的代价从而更关注少数类的学习。配合AUC评估基本能达到和小数据集过采样接近的效果而且不容易过拟合。如果你有足够多的数据上万条SMOTE的优势会逐渐显现但几千条以下还是保守一点比较好。6.3 单模型之外软投票集成的实际增量最后能不能再往上提一点性能我在单个模型达到瓶颈之后做了软投票集成Soft Voting。原理很简单把RF、SVM、KNN、神经网络四个分类器预测出的阳性概率取平均作为集成模型的输出。实际增量有多少拿ERα活性分类来说最好的单模型AUC是0.90软投票集成之后到了0.91到0.92之间。增幅不算大但胜在稳定而且预测概率更平滑不太容易出现某个模型在某个区域集体犯错的情况。更重要的是集成之后模型对随机种子的敏感度降低了这在项目交付时是非常重要的一点——单次结果波动大合作方会质疑你的模型稳定性。回归任务上的集成也类似把随机森林回归和SVR的预测值做加权平均通常能小幅提升R²并降低RMSE。加权系数可以用验证集上的表现来确定比如RF:RVR 0.6:0.4这样。6.4 实战中保留的其他经验细节再补充几个写代码时容易忽略的小细节。特征重要性分析随机森林模型可以直接输出feature_importances_。在ERα拮抗剂模型里我发现重要性排名靠前的往往是和芳环、疏水中心、氢键受体相关的指纹位这和ERα配体的已知构效关系高度吻合——配体结合口袋里有几个关键的极性残基所以有氢键受体的分子更容易有拮抗活性。这种可解释性分析不仅在论文里有用跟药化同事讨论后续优化方向时也非常有价值。PCA载荷分析也有类似作用。看前几个主成分里哪些原始特征权重最大可以对模型的特征意义有个整体感知。虽然不能直接说某个指纹位代表一个明确的官能团但结合Morgan指纹的半径范围去回溯子结构是可以验证模型的合理性。关于RDKit版本和随机种子我必须多说一句这个项目的可复现性比想象中重要。RDKit不同的版本生成的指纹在极少情况下会有细微差异scikit-learn/SVM的SVC在多线程情况下也可能有微小波动。所以实验开始前就把random_state42固定好并在记录实验结果的表格里写上环境版本号这是在给自己省未来的麻烦。用我自己的体会来收个尾跑通这个项目的技术栈并不复杂随便找个教程都能把代码跑起来。但这个项目的真正价值在于理解每个决策背后的权衡比如什么时候该用PCA什么时候该相信模型的AUC什么时候该对一个好得不真实的结果保持警惕。这些经验没有现成代码能给你只能在反复试错里积累。我希望这篇文章至少能让你在开始自己的药物数据预测项目时少走几步我当时绕过的弯路。