简介这份资源围绕除风清脾汤治疗血吸虫病的机制研究展开面向具备生物信息学与机器学习基础、从事中医药现代化研究的科研人员。内容整合网络药理学、机器学习、分子对接与分子动力学模拟完整呈现从TCMSP与UniProt获取成分及靶点、构建草药-靶点网络、Venn图筛选共同靶点到PPI分析、GO与KEGG富集、LASSO与随机森林及SVM-RFE筛选关键靶点再到汉黄芩素、山奈酚、木犀草素、槲皮素与靶点对接验证的全流程并附详细可运行代码与解释。资源包为1个PDF文件约808KB便于集中阅读与对照复现。目前已有90人学习下载。读者可借此掌握多成分-多靶点-多通路研究思路理解山奈酚与TP53稳定结合等关键发现并将代码与方法迁移至类似中药复方机制研究作为可参考的分析模板。1. 从一张药方到一套计算流程除风清脾汤治血吸虫病到底怎么研究除风清脾汤Chufeng Qingpi decoctionCQD是一张传统方剂而血吸虫病是由血吸虫寄生于门静脉系统引起的寄生虫病病理核心在于虫卵沉积诱发的肝脏肉芽肿与纤维化。把这两者放进同一句话里很多人第一反应是「中药复方成分那么杂靶点那么多怎么证明它有用」。这正是网络药理学、机器学习、分子对接和分子动力学模拟这套组合拳要回答的问题不是去证明某一味药杀死了成虫而是回答 CQD 里的哪些成分、通过哪些人体靶点、以什么结合方式干预了血吸虫病相关的炎症与纤维化通路。这套流程适合三类人做中药药理机制研究的研究生、想用 Python 把网络药理学跑通的计算生物学入门者、以及需要给复方机制找计算证据的临床科研人员。它不要求你会养虫、会做动物实验但要求你能装 Python、能看懂 SMILES 和 PDB 文件、能接受「计算结果只是假设最终还要实验验证」这个前提。下面按「数据从哪来 → 网络怎么建 → 机器学习怎么筛 → 对接和动力学怎么验」的顺序把每一步的命令、参数和翻车点讲清楚。2. 网络药理学打底CQD 成分与血吸虫病靶点怎么拿到手2.1 成分收集与 ADME 筛选的取舍CQD 的化学成分没有现成的单一权威清单常见做法是拆方检索把方中每味药分别到 TCMSP、HERB、SymMap 这类中药数据库里查再合并去重。检索时用拉丁名或中文名都试一遍因为不同库的命名规范不一致。拿到成分后第一件事是加 SMILES 和 PubChem CID没有 SMILES 的成分后面全部跑不动。ADME 筛选是网络药理学最容易被质疑的一步。经典标准是口服生物利用度 OB ≥ 30%、类药性 DL ≥ 0.18但这条线对苷类、多糖类成分极不友好很多已知有效成分会被直接筛掉。我的做法是双轨先用 OB≥30%、DL≥0.18 得到核心成分集再单独保留血吸虫病相关文献里明确报道过的成分两组合并。这样既符合审稿人对标准参数的期待又不至于把可能的关键成分误杀。import pandas as pd # 读取从 TCMSP 导出的成分表列名按实际导出调整 df pd.read_csv(cqd_ingredients_raw.csv) # 第一步去重同一成分可能来自多味药 df df.drop_duplicates(subset[MOL_ID]) # 第二步ADME 硬筛OB 和 DL 列必须是数值型 df[OB] pd.to_numeric(df[OB], errorscoerce) df[DL] pd.to_numeric(df[DL], errorscoerce) core df[(df[OB] 30) (df[DL] 0.18)].copy() # 第三步文献补充成分手动维护一个 MOL_ID 白名单 literature_ids [MOL000098, MOL000422] # 示例按实际文献替换 extra df[df[MOL_ID].isin(literature_ids)] final pd.concat([core, extra]).drop_duplicates(subset[MOL_ID]) final.to_csv(cqd_ingredients_final.csv, indexFalse) print(f核心成分 {len(core)} 个合并后 {len(final)} 个)这段代码的逻辑是「先标准筛、再人工补」errorscoerce是为了防止导出表里混入「-」或空字符串导致比较报错。参数上 OB 和 DL 的阈值可以调但一旦调了就要在方法学里写明理由否则审稿人会追问为什么是 25 不是 30。跑完这一步你会得到一个几十到一两百个成分的表这是后面所有网络的起点。2.2 靶点预测与血吸虫病靶点集的合并成分靶点预测常用 SwissTargetPrediction、TargetNet、PharmMapper。SwissTargetPrediction 一次最多提交几十个 SMILES批量跑要分批。血吸虫病靶点则从 GeneCards、OMIM、DisGeNET 搜「schistosomiasis」「liver fibrosis」「granuloma」等关键词注意 GeneCards 的 relevance score 要设个下限不然会混进大量弱相关基因。两边靶点拿到后必须做基因名标准化这是血泪经验最多的地方。不同库有的用 HGNC 官方 symbol有的用别名有的还是旧名。不统一的话后面取交集会凭空少掉一半。用 MyGene.info 或 org.Hs.eg.db 做映射把别名、旧 symbol 全部归一到当前官方 symbol。import mygene mg mygene.MyGeneInfo() def normalize_genes(gene_list): # 批量查询fields 指定返回官方 symbol res mg.querymany(gene_list, scopesalias,symbol,entrezgene, fieldssymbol, specieshuman) mapping {} for item in res: if symbol in item: mapping[item[query]] item[symbol] return mapping cqd_targets normalize_genes(cqd_gene_list) disease_targets normalize_genes(disease_gene_list) # 取交集得到潜在作用靶点 common set(cqd_targets.values()) set(disease_targets.values()) print(f交集靶点 {len(common)} 个)scopes里把 alias 和 symbol 都放进去是因为你手上的列表往往混着两种写法。specieshuman必须写否则会返回小鼠、大鼠的同源基因物种一乱整个网络就没意义了。交集靶点数量通常在几十到两百之间太少说明筛选太严太多说明疾病靶点没控好。2.3 构建 PPI 网络并找出真正该关注的核心靶点把交集靶点丢进 STRING物种选 Homo sapiensminimum interaction score 一般设 0.4medium confidence。下载 TSV 后在 Cytoscape 里做拓扑分析或者直接用 Python 的 networkx 算 degree、betweenness、closeness。核心靶点不是 degree 最高的那几个就完事通常要 degree 和 betweenness 双高才说明它既是枢纽又在信息流上关键。import networkx as nx G nx.from_pandas_edgelist(edges, sourcenode1, targetnode2) deg dict(G.degree()) btw nx.betweenness_centrality(G) # 按 degree 排序取前 15 作为候选核心 top_deg sorted(deg.items(), keylambda x: x[1], reverseTrue)[:15] for node, d in top_deg: print(node, d, round(btw[node], 4))betweenness 计算复杂度高节点上百个时用k100做近似采样不然会跑很久。这一步的输出就是后面分子对接的受体清单也是机器学习分类任务里的正样本来源之一。3. 机器学习上场把「可能有效」筛成「优先验证」3.1 为什么网络药理学之后还要加机器学习纯网络药理学的结论高度依赖数据库覆盖度和阈值设定换个库结果就变。机器学习在这里的作用不是替代网络分析而是给成分-靶点关系加一层基于已知活性数据的判别。常见做法是构建二分类模型正样本是已知对血吸虫或肝纤维化有活性的化合物负样本是确认无活性的化合物特征用分子指纹或描述符模型用随机森林、XGBoost 或 SVM。这一步的价值在于它能把网络药理学给出的几十个候选成分按预测概率重新排序让你优先去对接和做实验的那几个更有依据。但要注意正负样本的界定必须干净如果负样本里混进了没测过的化合物模型学到的就是噪声。3.2 用 RDKit 生成分子指纹并训练分类模型from rdkit import Chem from rdkit.Chem import AllChem import numpy as np from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score def smiles_to_fp(smi, n_bits2048): mol Chem.MolFromSmiles(smi) if mol is None: return None # Morgan 指纹半径 2 是常用默认值 fp AllChem.GetMorganFingerprintAsBitVect(mol, radius2, nBitsn_bits) return np.array(fp) # 构建特征矩阵跳过无法解析的 SMILES X, y [], [] for smi, label in zip(smiles_list, labels): fp smiles_to_fp(smi) if fp is not None: X.append(fp) y.append(label) X np.array(X) y np.array(y) clf RandomForestClassifier(n_estimators500, max_depthNone, class_weightbalanced, random_state42) scores cross_val_score(clf, X, y, cv5, scoringroc_auc) print(AUC:, scores.mean(), scores.std())radius2对应 ECFP4是活性预测里最常用的指纹。class_weightbalanced很重要因为活性化合物通常远少于非活性不加权模型会倾向于全预测负类AUC 看着还行但实际没用。n_estimators500是精度和速度的折中数据量小于一千时再往上加收益很小。交叉验证用 5 折如果 AUC 低于 0.7先别急着调参回去检查正负样本是不是有重叠或标签错误。3.3 特征重要性与模型可解释性随机森林能输出 feature importance但 Morgan 指纹的每一位对应什么子结构并不直观。想解释「模型到底看中了哪些结构」可以换成基于描述符的模型或者用 SHAP 值分析。实操里更省事的做法是把预测概率最高的前 20 个成分单独拿出来看它们有没有共同骨架再回到文献里查这类骨架有没有抗寄生虫或抗纤维化报道。这比强行解释每一位指纹更符合实际研究逻辑。import shap explainer shap.TreeExplainer(clf) shap_values explainer.shap_values(X[:200]) # 抽样解释全量太慢 # 对正类活性的贡献 shap.summary_plot(shap_values[1], X[:200])SHAP 在指纹这种高维稀疏特征上跑得慢抽样两百个样本足够看趋势。如果发现模型主要靠某几个位点判断而这些位点对应的子结构在正负样本里分布差异极大要警惕是不是数据泄漏——比如正样本全来自同一篇文献、同一类骨架。4. 分子对接把核心成分和核心靶点对上4.1 受体和配体的准备对接前受体要处理从 PDB 下载蛋白结构去水、去配体、加氢、加电荷。常用工具是 AutoDockTools 或 PyMOL 配合 prepare_receptor。配体就是前面筛出的成分用 RDKit 或 Open Babel 从 SMILES 生成 3D 构象并加氢。这一步的坑在于蛋白的质子化状态组氨酸、天冬氨酸、谷氨酸在不同 pH 下带电情况不同直接影响氢键和盐桥。没有实验结构时AlphaFold 预测结构可以用但要清楚它给的是静态构象侧链位置未必准。# 用 Open Babel 从 SMILES 生成 3D 配体加氢并做能量最小化 obabel -:CC(O)Oc1ccccc1C(O)O -O ligand.pdb --gen3d --minimize --ff MMFF94--gen3d生成三维坐标--minimize做初步优化--ff MMFF94指定力场。生成的构象质量直接影响对接结果如果配体是柔性大分子建议用 RDKit 生成多个构象再分别对接取最优打分。4.2 AutoDock Vina 对接与打分解读# 受体和配体都转成 pdbqt 后运行 Vina vina --receptor receptor.pdbqt --ligand ligand.pdbqt \ --center_x 10.0 --center_y 20.0 --center_z 15.0 \ --size_x 30 --size_y 30 --size_z 30 \ --exhaustiveness 32 --num_modes 9 --out out.pdbqt--center和--size定义对接盒子盒子要覆盖已知活性口袋或预测的结合位点太小会漏掉正确构象太大则计算量暴涨且打分区分度下降。--exhaustiveness 32比默认的 8 更充分适合最终验证阶段早期大批量筛选可以用 8 提速。打分单位是 kcal/mol负值越大结合越强但 Vina 打分和真实亲和力只是弱相关一般把 -7 kcal/mol 作为「值得进一步看」的参考线不要当成硬标准。对接结果要用 PyMOL 或 Discovery Studio 看相互作用氢键、π-π 堆积、疏水接触。如果核心成分和核心靶点的结合模式里出现了与已知抑制剂相似的关键残基相互作用这个结果的说服力会强很多。5. 分子动力学模拟对接结果的稳定性验证与常见翻车点5.1 模拟体系搭建与力场选择对接给的是静态快照分子动力学MD模拟看的是这个复合物在溶剂里跑一段时间后还稳不稳。常用 GROMACS 或 AMBER。蛋白力场选 AMBER99SB-ILDN 或 CHARMM36小分子配体力场用 GAFF2 配合 AM1-BCC 电荷这是目前复方成分模拟最通用的组合。# 用 GROMACS 做能量最小化 gmx grompp -f minim.mdp -c complex.gro -p topol.top -o em.tpr gmx mdrun -v -deffnm em # 平衡阶段NVT 和 NPT 各跑 100-500 ps gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx mdrun -v -deffnm nvtminim.mdp里最关键的是emtol和nsteps通常设 1000 kJ/mol/nm 和 50000 步。NVT 阶段用 V-rescale 控温 300 KNPT 阶段用 Parrinello-Rahman 控压 1 bar。盒子边界要保证蛋白距离盒子边缘至少 1.0 nm否则周期性镜像会让蛋白和自己打架。溶剂用 TIP3P 水加 0.15 mol/L NaCl 中和电荷并模拟生理离子强度。5.2 轨迹分析RMSD、RMSF 和结合自由能跑完 100 ns 生产模拟后先看 RMSD 判断体系是否平衡。蛋白骨架 RMSD 在 20 ns 后稳定在 0.2-0.3 nm 通常说明体系收敛。配体 RMSD 如果持续漂移说明结合不稳定对接结果可能有问题。# 计算蛋白骨架 RMSD gmx rms -s md.tpr -f md.xtc -o rmsd.xvg -tu ns # 计算每残基 RMSF gmx rmsf -s md.tpr -f md.xtc -o rmsf.xvg -res # MM-PBSA 结合自由能 gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx -cg 1 13 -ct md.xtcRMSF 高的残基是柔性区域如果结合口袋附近的残基 RMSF 在结合后明显降低说明配体让口袋变稳定了这是支持结合的证据。MM-PBSA 把结合自由能拆成范德华、静电、极性溶剂化和非极性溶剂化四项重点看范德华和静电的贡献如果这两项是主要驱动力和对接打分里的疏水、氢键结论能对上整个证据链就自洽了。5.3 避坑与排查分子模拟里最容易翻车的五件事现象一能量最小化不收敛报错提示力过大。原因通常是初始结构有原子重叠或者配体拓扑文件里的电荷、键参数不对。解决是先单独对配体做真空最小化再放回复合物检查 GAFF2 分配的原子类型是否合理必要时手动改。现象二NPT 平衡时盒子体积剧烈波动。多半是压力耦合参数设错或者体系里有真空区。检查tau_p是否在 2-5 ps 之间compressibility是否设成 4.5e-5。如果盒子是刚性的确认pcoupl没写成 no。现象三RMSD 一直不收敛跑 100 ns 还在涨。可能是蛋白有大的构象变化也可能是模拟时间不够。先延长到 200 ns 看趋势如果还在漂检查是不是 N 端或 C 端游离片段在乱动必要时对末端做位置限制。现象四配体跑出结合口袋。对接打分高不代表结合稳定。如果配体在 10 ns 内就离开口袋要么对接构象本身不合理要么力场参数让配体被溶剂拉走。回看对接时的关键相互作用有没有在模拟里维持没有的话这个靶点-成分对就要降级。现象五MM-PBSA 算出来结合自由能是正的。先别怀疑结论检查index.ndx里受体和配体的组号有没有选错这是最高频的低级错误。其次看熵项有没有算不算熵时自由能偏负是正常的但如果是正值多半是静电项被溶剂化抵消过头检查介电常数设置。6. 把整条链路串起来从结果到可验证假设的最后一公里跑完上面所有步骤你手上会有一张成分-靶点-通路网络、一个机器学习排序、若干对接构象和几条 MD 轨迹。真正决定这套研究值不值得做的是能不能从中提炼出两三个可实验验证的假设。我的习惯是做一个交叉优先级表成分在机器学习里预测概率高、在对接里打分低结合强、在 MD 里 RMSD 稳定、且靶点在 PPI 网络里 degree 和 betweenness 双高——同时满足这四条的组合优先送去做细胞或动物实验。证据层关注指标建议阈值不满足时的处理网络药理学靶点 degree / betweenness双高前 15降为背景靶点机器学习预测概率 0.7回查指纹特征分子对接Vina 打分 -7 kcal/mol换构象或换口袋MD 模拟配体 RMSD稳定在 0.15 nm 内延长模拟或弃用MM-PBSA范德华静电贡献为主要驱动力检查组号与熵项还有一个容易被忽略的技巧把 MD 轨迹里配体和靶点形成氢键的占有率算出来。占有率高于 60% 的氢键比对接图里画出来的单帧氢键可信得多。命令上可以用gmx hbond配合-hbm和-hbmc调距离和角度阈值默认 0.35 nm 和 30 度对大多数体系够用如果配体是柔性长链把角度放宽到 40 度再看。# 统计配体-蛋白氢键占有率 gmx hbond -s md.tpr -f md.xtc -n index.ndx -num hbnum.xvg # index.ndx 里需要提前把配体和蛋白定义成两个组最后说个我自己的教训早期做这类研究时我总想把所有成分、所有靶点、所有通路都塞进一张大图里结果图很漂亮但没人知道该验证什么。后来改成「每个结论都必须能落到一个具体的成分-靶点对和一个具体的实验读段上」文章的说服力和自己的实验效率都上来了。计算只是把搜索空间缩小真正拍板还得靠湿实验。希望帮到你。本文还有配套的精品资源点击获取 网络药理学+机器学习+分子对接与动力学：复方干预血吸虫病研究全流程
简介这份资源围绕除风清脾汤治疗血吸虫病的机制研究展开面向具备生物信息学与机器学习基础、从事中医药现代化研究的科研人员。内容整合网络药理学、机器学习、分子对接与分子动力学模拟完整呈现从TCMSP与UniProt获取成分及靶点、构建草药-靶点网络、Venn图筛选共同靶点到PPI分析、GO与KEGG富集、LASSO与随机森林及SVM-RFE筛选关键靶点再到汉黄芩素、山奈酚、木犀草素、槲皮素与靶点对接验证的全流程并附详细可运行代码与解释。资源包为1个PDF文件约808KB便于集中阅读与对照复现。目前已有90人学习下载。读者可借此掌握多成分-多靶点-多通路研究思路理解山奈酚与TP53稳定结合等关键发现并将代码与方法迁移至类似中药复方机制研究作为可参考的分析模板。1. 从一张药方到一套计算流程除风清脾汤治血吸虫病到底怎么研究除风清脾汤Chufeng Qingpi decoctionCQD是一张传统方剂而血吸虫病是由血吸虫寄生于门静脉系统引起的寄生虫病病理核心在于虫卵沉积诱发的肝脏肉芽肿与纤维化。把这两者放进同一句话里很多人第一反应是「中药复方成分那么杂靶点那么多怎么证明它有用」。这正是网络药理学、机器学习、分子对接和分子动力学模拟这套组合拳要回答的问题不是去证明某一味药杀死了成虫而是回答 CQD 里的哪些成分、通过哪些人体靶点、以什么结合方式干预了血吸虫病相关的炎症与纤维化通路。这套流程适合三类人做中药药理机制研究的研究生、想用 Python 把网络药理学跑通的计算生物学入门者、以及需要给复方机制找计算证据的临床科研人员。它不要求你会养虫、会做动物实验但要求你能装 Python、能看懂 SMILES 和 PDB 文件、能接受「计算结果只是假设最终还要实验验证」这个前提。下面按「数据从哪来 → 网络怎么建 → 机器学习怎么筛 → 对接和动力学怎么验」的顺序把每一步的命令、参数和翻车点讲清楚。2. 网络药理学打底CQD 成分与血吸虫病靶点怎么拿到手2.1 成分收集与 ADME 筛选的取舍CQD 的化学成分没有现成的单一权威清单常见做法是拆方检索把方中每味药分别到 TCMSP、HERB、SymMap 这类中药数据库里查再合并去重。检索时用拉丁名或中文名都试一遍因为不同库的命名规范不一致。拿到成分后第一件事是加 SMILES 和 PubChem CID没有 SMILES 的成分后面全部跑不动。ADME 筛选是网络药理学最容易被质疑的一步。经典标准是口服生物利用度 OB ≥ 30%、类药性 DL ≥ 0.18但这条线对苷类、多糖类成分极不友好很多已知有效成分会被直接筛掉。我的做法是双轨先用 OB≥30%、DL≥0.18 得到核心成分集再单独保留血吸虫病相关文献里明确报道过的成分两组合并。这样既符合审稿人对标准参数的期待又不至于把可能的关键成分误杀。import pandas as pd # 读取从 TCMSP 导出的成分表列名按实际导出调整 df pd.read_csv(cqd_ingredients_raw.csv) # 第一步去重同一成分可能来自多味药 df df.drop_duplicates(subset[MOL_ID]) # 第二步ADME 硬筛OB 和 DL 列必须是数值型 df[OB] pd.to_numeric(df[OB], errorscoerce) df[DL] pd.to_numeric(df[DL], errorscoerce) core df[(df[OB] 30) (df[DL] 0.18)].copy() # 第三步文献补充成分手动维护一个 MOL_ID 白名单 literature_ids [MOL000098, MOL000422] # 示例按实际文献替换 extra df[df[MOL_ID].isin(literature_ids)] final pd.concat([core, extra]).drop_duplicates(subset[MOL_ID]) final.to_csv(cqd_ingredients_final.csv, indexFalse) print(f核心成分 {len(core)} 个合并后 {len(final)} 个)这段代码的逻辑是「先标准筛、再人工补」errorscoerce是为了防止导出表里混入「-」或空字符串导致比较报错。参数上 OB 和 DL 的阈值可以调但一旦调了就要在方法学里写明理由否则审稿人会追问为什么是 25 不是 30。跑完这一步你会得到一个几十到一两百个成分的表这是后面所有网络的起点。2.2 靶点预测与血吸虫病靶点集的合并成分靶点预测常用 SwissTargetPrediction、TargetNet、PharmMapper。SwissTargetPrediction 一次最多提交几十个 SMILES批量跑要分批。血吸虫病靶点则从 GeneCards、OMIM、DisGeNET 搜「schistosomiasis」「liver fibrosis」「granuloma」等关键词注意 GeneCards 的 relevance score 要设个下限不然会混进大量弱相关基因。两边靶点拿到后必须做基因名标准化这是血泪经验最多的地方。不同库有的用 HGNC 官方 symbol有的用别名有的还是旧名。不统一的话后面取交集会凭空少掉一半。用 MyGene.info 或 org.Hs.eg.db 做映射把别名、旧 symbol 全部归一到当前官方 symbol。import mygene mg mygene.MyGeneInfo() def normalize_genes(gene_list): # 批量查询fields 指定返回官方 symbol res mg.querymany(gene_list, scopesalias,symbol,entrezgene, fieldssymbol, specieshuman) mapping {} for item in res: if symbol in item: mapping[item[query]] item[symbol] return mapping cqd_targets normalize_genes(cqd_gene_list) disease_targets normalize_genes(disease_gene_list) # 取交集得到潜在作用靶点 common set(cqd_targets.values()) set(disease_targets.values()) print(f交集靶点 {len(common)} 个)scopes里把 alias 和 symbol 都放进去是因为你手上的列表往往混着两种写法。specieshuman必须写否则会返回小鼠、大鼠的同源基因物种一乱整个网络就没意义了。交集靶点数量通常在几十到两百之间太少说明筛选太严太多说明疾病靶点没控好。2.3 构建 PPI 网络并找出真正该关注的核心靶点把交集靶点丢进 STRING物种选 Homo sapiensminimum interaction score 一般设 0.4medium confidence。下载 TSV 后在 Cytoscape 里做拓扑分析或者直接用 Python 的 networkx 算 degree、betweenness、closeness。核心靶点不是 degree 最高的那几个就完事通常要 degree 和 betweenness 双高才说明它既是枢纽又在信息流上关键。import networkx as nx G nx.from_pandas_edgelist(edges, sourcenode1, targetnode2) deg dict(G.degree()) btw nx.betweenness_centrality(G) # 按 degree 排序取前 15 作为候选核心 top_deg sorted(deg.items(), keylambda x: x[1], reverseTrue)[:15] for node, d in top_deg: print(node, d, round(btw[node], 4))betweenness 计算复杂度高节点上百个时用k100做近似采样不然会跑很久。这一步的输出就是后面分子对接的受体清单也是机器学习分类任务里的正样本来源之一。3. 机器学习上场把「可能有效」筛成「优先验证」3.1 为什么网络药理学之后还要加机器学习纯网络药理学的结论高度依赖数据库覆盖度和阈值设定换个库结果就变。机器学习在这里的作用不是替代网络分析而是给成分-靶点关系加一层基于已知活性数据的判别。常见做法是构建二分类模型正样本是已知对血吸虫或肝纤维化有活性的化合物负样本是确认无活性的化合物特征用分子指纹或描述符模型用随机森林、XGBoost 或 SVM。这一步的价值在于它能把网络药理学给出的几十个候选成分按预测概率重新排序让你优先去对接和做实验的那几个更有依据。但要注意正负样本的界定必须干净如果负样本里混进了没测过的化合物模型学到的就是噪声。3.2 用 RDKit 生成分子指纹并训练分类模型from rdkit import Chem from rdkit.Chem import AllChem import numpy as np from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score def smiles_to_fp(smi, n_bits2048): mol Chem.MolFromSmiles(smi) if mol is None: return None # Morgan 指纹半径 2 是常用默认值 fp AllChem.GetMorganFingerprintAsBitVect(mol, radius2, nBitsn_bits) return np.array(fp) # 构建特征矩阵跳过无法解析的 SMILES X, y [], [] for smi, label in zip(smiles_list, labels): fp smiles_to_fp(smi) if fp is not None: X.append(fp) y.append(label) X np.array(X) y np.array(y) clf RandomForestClassifier(n_estimators500, max_depthNone, class_weightbalanced, random_state42) scores cross_val_score(clf, X, y, cv5, scoringroc_auc) print(AUC:, scores.mean(), scores.std())radius2对应 ECFP4是活性预测里最常用的指纹。class_weightbalanced很重要因为活性化合物通常远少于非活性不加权模型会倾向于全预测负类AUC 看着还行但实际没用。n_estimators500是精度和速度的折中数据量小于一千时再往上加收益很小。交叉验证用 5 折如果 AUC 低于 0.7先别急着调参回去检查正负样本是不是有重叠或标签错误。3.3 特征重要性与模型可解释性随机森林能输出 feature importance但 Morgan 指纹的每一位对应什么子结构并不直观。想解释「模型到底看中了哪些结构」可以换成基于描述符的模型或者用 SHAP 值分析。实操里更省事的做法是把预测概率最高的前 20 个成分单独拿出来看它们有没有共同骨架再回到文献里查这类骨架有没有抗寄生虫或抗纤维化报道。这比强行解释每一位指纹更符合实际研究逻辑。import shap explainer shap.TreeExplainer(clf) shap_values explainer.shap_values(X[:200]) # 抽样解释全量太慢 # 对正类活性的贡献 shap.summary_plot(shap_values[1], X[:200])SHAP 在指纹这种高维稀疏特征上跑得慢抽样两百个样本足够看趋势。如果发现模型主要靠某几个位点判断而这些位点对应的子结构在正负样本里分布差异极大要警惕是不是数据泄漏——比如正样本全来自同一篇文献、同一类骨架。4. 分子对接把核心成分和核心靶点对上4.1 受体和配体的准备对接前受体要处理从 PDB 下载蛋白结构去水、去配体、加氢、加电荷。常用工具是 AutoDockTools 或 PyMOL 配合 prepare_receptor。配体就是前面筛出的成分用 RDKit 或 Open Babel 从 SMILES 生成 3D 构象并加氢。这一步的坑在于蛋白的质子化状态组氨酸、天冬氨酸、谷氨酸在不同 pH 下带电情况不同直接影响氢键和盐桥。没有实验结构时AlphaFold 预测结构可以用但要清楚它给的是静态构象侧链位置未必准。# 用 Open Babel 从 SMILES 生成 3D 配体加氢并做能量最小化 obabel -:CC(O)Oc1ccccc1C(O)O -O ligand.pdb --gen3d --minimize --ff MMFF94--gen3d生成三维坐标--minimize做初步优化--ff MMFF94指定力场。生成的构象质量直接影响对接结果如果配体是柔性大分子建议用 RDKit 生成多个构象再分别对接取最优打分。4.2 AutoDock Vina 对接与打分解读# 受体和配体都转成 pdbqt 后运行 Vina vina --receptor receptor.pdbqt --ligand ligand.pdbqt \ --center_x 10.0 --center_y 20.0 --center_z 15.0 \ --size_x 30 --size_y 30 --size_z 30 \ --exhaustiveness 32 --num_modes 9 --out out.pdbqt--center和--size定义对接盒子盒子要覆盖已知活性口袋或预测的结合位点太小会漏掉正确构象太大则计算量暴涨且打分区分度下降。--exhaustiveness 32比默认的 8 更充分适合最终验证阶段早期大批量筛选可以用 8 提速。打分单位是 kcal/mol负值越大结合越强但 Vina 打分和真实亲和力只是弱相关一般把 -7 kcal/mol 作为「值得进一步看」的参考线不要当成硬标准。对接结果要用 PyMOL 或 Discovery Studio 看相互作用氢键、π-π 堆积、疏水接触。如果核心成分和核心靶点的结合模式里出现了与已知抑制剂相似的关键残基相互作用这个结果的说服力会强很多。5. 分子动力学模拟对接结果的稳定性验证与常见翻车点5.1 模拟体系搭建与力场选择对接给的是静态快照分子动力学MD模拟看的是这个复合物在溶剂里跑一段时间后还稳不稳。常用 GROMACS 或 AMBER。蛋白力场选 AMBER99SB-ILDN 或 CHARMM36小分子配体力场用 GAFF2 配合 AM1-BCC 电荷这是目前复方成分模拟最通用的组合。# 用 GROMACS 做能量最小化 gmx grompp -f minim.mdp -c complex.gro -p topol.top -o em.tpr gmx mdrun -v -deffnm em # 平衡阶段NVT 和 NPT 各跑 100-500 ps gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx mdrun -v -deffnm nvtminim.mdp里最关键的是emtol和nsteps通常设 1000 kJ/mol/nm 和 50000 步。NVT 阶段用 V-rescale 控温 300 KNPT 阶段用 Parrinello-Rahman 控压 1 bar。盒子边界要保证蛋白距离盒子边缘至少 1.0 nm否则周期性镜像会让蛋白和自己打架。溶剂用 TIP3P 水加 0.15 mol/L NaCl 中和电荷并模拟生理离子强度。5.2 轨迹分析RMSD、RMSF 和结合自由能跑完 100 ns 生产模拟后先看 RMSD 判断体系是否平衡。蛋白骨架 RMSD 在 20 ns 后稳定在 0.2-0.3 nm 通常说明体系收敛。配体 RMSD 如果持续漂移说明结合不稳定对接结果可能有问题。# 计算蛋白骨架 RMSD gmx rms -s md.tpr -f md.xtc -o rmsd.xvg -tu ns # 计算每残基 RMSF gmx rmsf -s md.tpr -f md.xtc -o rmsf.xvg -res # MM-PBSA 结合自由能 gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx -cg 1 13 -ct md.xtcRMSF 高的残基是柔性区域如果结合口袋附近的残基 RMSF 在结合后明显降低说明配体让口袋变稳定了这是支持结合的证据。MM-PBSA 把结合自由能拆成范德华、静电、极性溶剂化和非极性溶剂化四项重点看范德华和静电的贡献如果这两项是主要驱动力和对接打分里的疏水、氢键结论能对上整个证据链就自洽了。5.3 避坑与排查分子模拟里最容易翻车的五件事现象一能量最小化不收敛报错提示力过大。原因通常是初始结构有原子重叠或者配体拓扑文件里的电荷、键参数不对。解决是先单独对配体做真空最小化再放回复合物检查 GAFF2 分配的原子类型是否合理必要时手动改。现象二NPT 平衡时盒子体积剧烈波动。多半是压力耦合参数设错或者体系里有真空区。检查tau_p是否在 2-5 ps 之间compressibility是否设成 4.5e-5。如果盒子是刚性的确认pcoupl没写成 no。现象三RMSD 一直不收敛跑 100 ns 还在涨。可能是蛋白有大的构象变化也可能是模拟时间不够。先延长到 200 ns 看趋势如果还在漂检查是不是 N 端或 C 端游离片段在乱动必要时对末端做位置限制。现象四配体跑出结合口袋。对接打分高不代表结合稳定。如果配体在 10 ns 内就离开口袋要么对接构象本身不合理要么力场参数让配体被溶剂拉走。回看对接时的关键相互作用有没有在模拟里维持没有的话这个靶点-成分对就要降级。现象五MM-PBSA 算出来结合自由能是正的。先别怀疑结论检查index.ndx里受体和配体的组号有没有选错这是最高频的低级错误。其次看熵项有没有算不算熵时自由能偏负是正常的但如果是正值多半是静电项被溶剂化抵消过头检查介电常数设置。6. 把整条链路串起来从结果到可验证假设的最后一公里跑完上面所有步骤你手上会有一张成分-靶点-通路网络、一个机器学习排序、若干对接构象和几条 MD 轨迹。真正决定这套研究值不值得做的是能不能从中提炼出两三个可实验验证的假设。我的习惯是做一个交叉优先级表成分在机器学习里预测概率高、在对接里打分低结合强、在 MD 里 RMSD 稳定、且靶点在 PPI 网络里 degree 和 betweenness 双高——同时满足这四条的组合优先送去做细胞或动物实验。证据层关注指标建议阈值不满足时的处理网络药理学靶点 degree / betweenness双高前 15降为背景靶点机器学习预测概率 0.7回查指纹特征分子对接Vina 打分 -7 kcal/mol换构象或换口袋MD 模拟配体 RMSD稳定在 0.15 nm 内延长模拟或弃用MM-PBSA范德华静电贡献为主要驱动力检查组号与熵项还有一个容易被忽略的技巧把 MD 轨迹里配体和靶点形成氢键的占有率算出来。占有率高于 60% 的氢键比对接图里画出来的单帧氢键可信得多。命令上可以用gmx hbond配合-hbm和-hbmc调距离和角度阈值默认 0.35 nm 和 30 度对大多数体系够用如果配体是柔性长链把角度放宽到 40 度再看。# 统计配体-蛋白氢键占有率 gmx hbond -s md.tpr -f md.xtc -n index.ndx -num hbnum.xvg # index.ndx 里需要提前把配体和蛋白定义成两个组最后说个我自己的教训早期做这类研究时我总想把所有成分、所有靶点、所有通路都塞进一张大图里结果图很漂亮但没人知道该验证什么。后来改成「每个结论都必须能落到一个具体的成分-靶点对和一个具体的实验读段上」文章的说服力和自己的实验效率都上来了。计算只是把搜索空间缩小真正拍板还得靠湿实验。希望帮到你。本文还有配套的精品资源点击获取