简介这份资源面向具备生物信息学、网络药理学或中医药研究背景的科研人员以及希望掌握中药复方机制研究流程的学者提供除风清脾汤治疗血吸虫病机制研究的完整论文复现资料。内容涵盖TCMSP与UniProt靶点识别、草药-靶点网络构建、Venn共同靶点筛选、PPI网络分析、GO与KEGG富集以及LASSO、随机森林、SVM-RFE等机器学习算法筛选关键靶点并通过分子对接与动力学模拟验证汉黄芩素、山奈酚等成分与TP53、TNF、IL6的相互作用。资源包为1个PDF文件大小约808KB内含可运行代码及逐段解释便于读者对照复现研究流程、理解多成分-多靶点-多通路机制并作为类似中药复方研究的参考模板。目前已有141人学习适合需要快速上手网络药理学与机器学习联合分析的研究者。1. 从一篇血吸虫病机制论文说起这套复现资源到底能跑出什么血吸虫病这块很多人第一反应是吡喹酮但耐药性和再感染的问题一直没消停所以中药复方作为替代或辅助方案的研究这几年明显多了起来。除风清脾汤CQD就是其中一个被拿来研究的方向问题是中药复方成分复杂到底哪个成分、作用在哪个靶点、走哪条通路光靠体外实验根本说不清。这套复现资源干的事就是用网络药理学筛靶点、机器学习做特征筛选和分类验证、分子对接看结合能不能成立、分子动力学模拟再确认复合物稳不稳把一条从成分到机制的完整证据链用代码串起来。适合谁做中药网络药理学的研究生、想补分子模拟实操的从业者以及需要一套完整复现流程当模板的人。它不是一个能一键出结果的工具包而是一套需要你理解每一步在干什么、参数为什么这么设的流程代码。2. 网络药理学打底成分-靶点-通路三张表怎么建2.1 为什么先做网络药理学而不是直接上分子对接中药复方的研究有个绕不开的坎你不可能把整锅汤拿去做对接必须先知道里面有哪些成分、这些成分可能作用在哪些靶点上。网络药理学的价值就在这——它用数据库检索和预测算法把化学成分和蛋白靶点之间的可能关系先铺出来形成一个候选池。没有这一步你后面做分子对接就是盲选选哪个蛋白、选哪个配体全靠拍脑袋审稿人第一个问题就能把你问住。常见做法是先从 TCMSP、BATMAN-TCM 这类数据库拉 CQD 的化学成分按 OB口服生物利用度≥30% 和 DL类药性≥0.18 筛一轮得到活性成分再用 SwissTargetPrediction 或类似工具预测成分靶点同时从 GeneCards、OMIM、DisGeNET 拉血吸虫病的疾病靶点两边取交集得到成分-疾病共同靶点。这套资源里的代码就是把这个流程脚本化了省得你一个个手动查。2.2 成分筛选与靶点预测的代码实现import pandas as pd # 读取 TCMSP 导出的成分表列名按实际导出调整 herbs pd.read_csv(cqd_components.csv) # 按 OB 和 DL 双阈值筛选活性成分 active herbs[(herbs[OB] 30) (herbs[DL] 0.18)].copy() print(f筛选前成分数: {len(herbs)}, 筛选后: {len(active)}) # 去重同一成分可能来自多味药 active_unique active.drop_duplicates(subset[MOL_ID, Molecule_name]) active_unique.to_csv(active_components.csv, indexFalse) # 读取靶点预测结果SwissTargetPrediction 导出格式 targets pd.read_csv(predicted_targets.csv) # 只保留 probability 0 的预测结果 targets targets[targets[Probability] 0] targets.to_csv(filtered_targets.csv, indexFalse)这段代码的逻辑很直白第一步按 OB 和 DL 卡阈值这是网络药理学里最常规的筛选标准OB 太低意味着吃进去吸收不了DL 太低意味着不像药。第二步去重是因为复方里不同药材可能含相同成分不去重后面做网络图会重复计数。第三步过滤靶点预测结果probability 为 0 的基本可以当噪声扔掉。参数方面OB 和 DL 的阈值不是死的有研究用 OB≥20% 或 DL≥0.10取决于你的复方类型和审稿要求。我一般会跑两套阈值对比如果结果差异大说明你的成分对阈值敏感需要在文章里说明选择理由。靶点预测的 probability 阈值同理SwissTargetPrediction 默认给的是 0 到 1 的概率值卡 0 只是去掉完全没预测到的实际分析中很多人会卡 0.1 或更高。2.3 交集靶点与 PPI 网络的构建拿到成分靶点和疾病靶点后取交集是标准操作。交集靶点才是后续 PPI 网络、富集分析和分子对接的输入。# 成分靶点与疾病靶点取交集 component_targets set(targets[Target]) disease_targets set(pd.read_csv(disease_targets.csv)[Gene]) common_targets component_targets disease_targets print(f成分靶点数: {len(component_targets)}, 疾病靶点数: {len(disease_targets)}) print(f交集靶点数: {len(common_targets)}) # 输出交集靶点列表用于 STRING 数据库构建 PPI pd.Series(list(common_targets)).to_csv(common_targets.txt, indexFalse, headerFalse)交集靶点拿到后丢进 STRING 数据库做 PPI 网络导出 TSV 文件再用 Cytoscape 做可视化和拓扑分析。这里有个容易翻车的地方STRING 的物种设置一定要选对人源Homo sapiens因为你的靶点预测结果大概率是人源蛋白如果物种选错PPI 网络会少一大截节点。拓扑分析里 degree、betweenness、closeness 三个指标通常用来筛核心靶点degree 排名前 10 或前 15 的节点就是后续分子对接的首选对象。3. 机器学习做特征筛选与分类验证别把靶点当特征硬塞3.1 机器学习在这个流程里到底扮演什么角色网络药理学给你的是一个候选靶点池但池子里靶点太多哪些才是真正跟血吸虫病相关的核心靶点这时候机器学习就派上用场了。常见做法是用交集靶点的表达数据比如 GEO 数据库里的血吸虫病相关数据集作为特征用 LASSO 回归或随机森林做特征筛选把冗余的、贡献小的靶点剔掉留下真正有区分度的特征。然后再用这些特征训练分类模型验证它们能不能把疾病样本和正常样本分开。这里要强调一点机器学习不是用来证明机制的它是用来做特征筛选和辅助验证的。你不可能用机器学习直接说CQD 通过靶点 X 治疗血吸虫病但你可以说通过 LASSO 筛选出 5 个核心靶点基于这 5 个靶点构建的随机森林模型在训练集和验证集上 AUC 均大于 0.85这是两个不同层面的结论。3.2 LASSO 特征筛选与随机森林分类的代码实现from sklearn.linear_model import LassoCV from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score, train_test_split from sklearn.metrics import roc_auc_score import numpy as np # X: 样本×靶点表达矩阵, y: 标签(1疾病, 0正常) X pd.read_csv(expression_matrix.csv, index_col0) y pd.read_csv(labels.csv, index_col0)[label] # LASSO 交叉验证筛选特征 lasso LassoCV(cv5, random_state42, max_iter10000) lasso.fit(X, y) # 保留系数非零的特征 selected X.columns[lasso.coef_ ! 0] print(fLASSO 筛选后保留特征数: {len(selected)}) print(f保留特征: {list(selected)}) # 用筛选后的特征训练随机森林 X_sel X[selected] X_train, X_test, y_train, y_test train_test_split( X_sel, y, test_size0.3, random_state42, stratifyy ) rf RandomForestClassifier(n_estimators500, random_state42) rf.fit(X_train, y_train) y_prob rf.predict_proba(X_test)[:, 1] print(f测试集 AUC: {roc_auc_score(y_test, y_prob):.3f}) # 交叉验证评估稳定性 cv_scores cross_val_score(rf, X_sel, y, cv5, scoringroc_auc) print(f5折交叉验证 AUC: {cv_scores.mean():.3f} ± {cv_scores.std():.3f})LASSO 那一步的核心参数是cv5表示 5 折交叉验证选正则化强度max_iter设大一点防止不收敛。系数非零的特征就是被保留的系数为零的被压缩掉了。随机森林这边n_estimators500是树的数量一般 300 到 1000 之间太少不稳定太多计算慢但收益递减。stratifyy保证训练集和测试集的类别比例一致样本不平衡时这个参数很关键。3.3 特征重要性排序与模型可解释性随机森林自带特征重要性可以直接拿来排序看哪些靶点对分类贡献最大。# 特征重要性排序 importance pd.Series(rf.feature_importances_, indexselected) importance importance.sort_values(ascendingFalse) print(importance.head(10)) # 输出用于后续分子对接的核心靶点 core_targets importance.head(5).index.tolist() print(f核心靶点: {core_targets})这里选前 5 个还是前 10 个没有硬性标准。我一般会结合 PPI 网络的 degree 排名一起看两个方法都排在前面的靶点可信度更高。如果 LASSO 筛完只剩两三个特征那说明你的数据维度可能有问题要么样本量太少要么靶点表达矩阵噪声太大这时候硬跑模型没意义得回头检查数据质量。注意机器学习模型的 AUC 高不代表机制成立它只说明这些靶点在表达层面有区分度。审稿人如果问你怎么证明 CQD 作用于这些靶点你还是得靠分子对接和动力学模拟来回答。4. 分子对接与分子动力学模拟从静态结合到动态稳定性4.1 分子对接的输入准备与参数设置分子对接干的事是把配体CQD 的活性成分放进受体核心靶点蛋白的活性口袋里算一个结合能看它们能不能稳定结合。这套流程里受体蛋白结构从 PDB 数据库下载配体结构从 PubChem 下载或自己构建。常见做法是用 AutoDock Vina 做对接因为它是开源的、速度也还行。准备工作包括受体去水、加氢、算电荷配体加氢、设可旋转键。这些用 AutoDockTools 或 Meeko 都能做。# 受体准备去水、加氢用 AutoDockTools 命令行或 PyMOL # 配体准备用 Meeko 转换格式 mk_prepare_ligand.py -i ligand.sdf -o ligand.pdbqt # 运行 AutoDock Vina 对接 vina --receptor receptor.pdbqt --ligand ligand.pdbqt \ --center_x 10.0 --center_y 20.0 --center_z 30.0 \ --size_x 20 --size_y 20 --size_z 20 \ --exhaustiveness 32 --num_modes 9 \ --out docking_result.pdbqt--center_x/y/z是活性口袋的中心坐标这个不能瞎设得从文献里查或者用工具预测。--size_x/y/z是搜索盒子的大小20×20×20 埃是个常用起点如果口袋比较大可以适当放大。--exhaustiveness 32是搜索彻底程度默认是 8调到 32 会慢一些但结果更稳。--num_modes 9是输出多少个构象一般看前 3 个就够了。对接结果里结合能affinity越低越好通常小于 -5.0 kcal/mol 认为有结合可能小于 -7.0 认为结合较强。但别只看结合能还要看结合模式——配体是不是真的进了口袋有没有跟关键残基形成氢键这些用 PyMOL 或 Discovery Studio 可视化确认。4.2 分子动力学模拟的体系构建与参数对接是静态的分子动力学模拟是动态的。它把对接得到的复合物放进一个模拟盒子里加水、加离子然后在一定的温度和压力下跑一段时间看复合物稳不稳定。RMSD 是衡量稳定性的核心指标如果 RMSD 在跑的过程中一直漂说明复合物不稳定对接结果可能不可靠。# 用 GROMACS 做分子动力学模拟的典型流程 # 1. 生成拓扑文件 gmx pdb2gmx -f complex.pdb -o complex.gro -water tip3p -ff amber99sb-ildn # 2. 定义模拟盒子 gmx editconf -f complex.gro -o box.gro -c -d 1.0 -bt cubic # 3. 加水 gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top # 4. 加离子中和体系 gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr gmx genion -s ions.tpr -o ionized.gro -p topol.top -pname NA -nname CL -neutral # 5. 能量最小化 gmx grompp -f minim.mdp -c ionized.gro -p topol.top -o em.tpr gmx mdrun -v -deffnm em # 6. NVT 平衡100ps gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx mdrun -deffnm nvt # 7. NPT 平衡100ps gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -p topol.top -o npt.tpr gmx mdrun -deffnm npt # 8. 成品模拟100ns gmx grompp -f md.mdp -c npt.gro -p topol.top -o md.tpr gmx mdrun -deffnm md这套流程里pdb2gmx选力场是关键amber99sb-ildn 是蛋白模拟里常用的-water tip3p是水模型。editconf的-d 1.0表示蛋白距离盒子边缘至少 1.0 nm太近会导致周期性边界效应。genion加离子中和体系电荷这一步不能省带电体系跑模拟会出问题。NVT 和 NPT 平衡各 100ps 是常规操作成品模拟跑 100ns 是很多论文的标配但如果你体系大、资源有限50ns 也能发文章关键看 RMSD 有没有收敛。4.3 RMSD、RMSF 与结合自由能计算模拟跑完后用gmx rms算 RMSDgmx rmsf算 RMSF看复合物稳不稳定、哪些残基波动大。结合自由能用 MM-PBSA 或 MM-GBSA 算gmx_MMPBSA 是常用的工具。# 计算 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 # 用 gmx_MMPBSA 算结合自由能 gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx -cg 1 13 -ct md.xtcRMSD 图如果在前 10ns 有上升然后平台说明体系平衡了。如果一直上升说明模拟时间不够或者体系本身不稳定。RMSF 看的是每个残基的波动活性口袋区域的残基如果 RMSF 很低说明结合后刚性好这是好现象。MM-PBSA 算出来的结合自由能如果是负的跟对接结果一致那这条证据链就比较完整了。5. 避坑与排查这套流程里最容易翻车的五个地方5.1 成分筛选阈值设太死活性成分漏掉一半现象按 OB≥30%、DL≥0.18 筛完CQD 只剩十几个成分文献里报道的某些已知活性成分没进去。原因OB 和 DL 的预测模型对某些结构类型比如苷类、多糖本身就不准阈值卡太死会误杀。解决跑两套阈值对比一套严格OB≥30%、DL≥0.18一套宽松OB≥20%、DL≥0.10取并集做后续分析在文章里说明阈值选择依据。如果某个成分文献明确报道有活性但被筛掉了可以手动加回来但要在方法里写清楚。5.2 靶点预测结果不统一不同数据库差出一大截现象SwissTargetPrediction 预测的靶点和 TCMSP 自带的靶点列表对不上交集少得可怜。原因不同数据库的预测算法和覆盖范围不一样SwissTargetPrediction 基于相似性TCMSP 基于实验数据结果有差异是正常的。解决不要只用一个数据库的靶点至少用两个数据库取并集然后再跟疾病靶点取交集。如果并集后靶点太多再用 PPI 网络的 degree 排名筛一轮。文章里要写清楚用了哪些数据库、怎么合并的。5.3 分子对接盒子设错位置配体根本没进活性口袋现象对接结果结合能很低比如 -9.0但可视化一看配体在蛋白表面根本没进活性口袋。原因--center_x/y/z设错了搜索盒子没覆盖真正的活性位点。解决活性口袋坐标从文献查或者用 CASTp、fpocket 预测。对接完必须可视化检查别只看结合能。如果配体在表面结合能再低也没意义。5.4 分子动力学模拟体系带电跑一半崩了现象模拟跑到几纳秒就报错提示体系不中性或者能量异常。原因genion那一步没做或者离子数量不对体系带净电荷。解决gmx genion的-neutral参数会自动中和体系但前提是ions.mdp里设了正确的离子类型。跑之前用gmx check检查体系电荷。另外如果配体带电荷pdb2gmx那一步要选对质子化状态。5.5 机器学习模型过拟合交叉验证 AUC 高但测试集崩了现象训练集 AUC 0.99交叉验证 AUC 0.95但独立测试集 AUC 只有 0.6。原因样本量太少、特征太多模型记住了训练数据的噪声。解决先做特征筛选降维LASSO 或随机森林都可以。样本量小于 50 的时候别用深度学习老老实实用 LASSO 逻辑回归或随机森林。交叉验证要用分层抽样测试集要独立于训练集。如果 AUC 还是崩说明数据本身有问题不是调参能解决的。6. 把复现流程跑通之后我习惯用这套检查清单收尾跑通一遍不代表结果可信我自己的习惯是每次复现完按下面这套清单过一遍确认没有低级错误再往下写。检查项具体操作通过标准成分筛选跑两套阈值对比活性成分数差异不超过 30%靶点交集至少两个数据库取并集交集靶点数在合理范围通常 50-200PPI 网络STRING 物种设为人源网络节点数与交集靶点数一致机器学习交叉验证 独立测试集测试集 AUC 与交叉验证 AUC 差距小于 0.1分子对接可视化检查结合模式配体在活性口袋内与关键残基有相互作用动力学模拟RMSD 收敛后 50ns RMSD 波动小于 0.2 nm结合自由能MM-PBSA 与对接结果对比趋势一致结合能均为负这套清单里分子对接的可视化检查和动力学模拟的 RMSD 收敛是我最看重的两步。对接结合能再漂亮配体没进袋子就是白搭RMSD 不收敛后面算的自由能都不可信。机器学习那块如果测试集 AUC 和交叉验证差太多别硬着头皮往下写回头查数据泄露或者样本标签有没有搞错。还有一个进阶用法把分子动力学模拟的轨迹做聚类分析看配体在口袋里的优势构象然后拿优势构象再跑一次 MM-PBSA这样算出来的结合自由能比单帧更可靠。gmx cluster可以做这件事选-method gromos按 RMSD 聚类看最大的簇对应的构象。从那以后我每次跑完分子动力学都强制自己先看 RMSD 和 RMSF 再往下走不收敛就加时间重跑绝不硬算自由能。希望帮到你。本文还有配套的精品资源点击获取