简介这份PDF资源面向具备生物信息学、网络药理学或中医药研究背景的科研人员聚焦除风清脾汤治疗血吸虫病的机制复现。内容完整呈现从TCMSP与UniProt识别成分靶点、构建草药-靶点网络到Venn图取交集、PPI分析、GO与KEGG富集再到LASSO、随机森林、SVM-RFE筛选关键靶点最后以分子对接和动力学模拟验证汉黄芩素、山奈酚等成分与TP53、TNF、IL6结合稳定性的全流程。资源包仅1个PDF文件约808KB内含可运行Python代码及逐段解释便于对照复现。已有141人学习。读者可借此掌握多成分-多靶点-多通路研究范式理解机器学习与网络药理学结合筛选疾病靶点的思路并获取分子对接与动力学验证的代码模板和排错参考适合作为类似中药复方机制研究的复现范例。1. 除风清脾汤治血吸虫病一篇论文复现为什么值得动手做血吸虫病这个方向很多人第一反应是“离我太远”。但如果你正在做网络药理学、机器学习或者分子动力学模拟的课题这篇以除风清脾汤CQD为对象的复现工作恰好是一个把四套方法串成完整证据链的样本。它要回答的问题很具体一味由多味中药组成的复方究竟通过哪些成分、哪些靶点、哪条通路对血吸虫病产生作用。单靠网络药理学只能给出“可能相关”的候选清单单靠分子对接只能证明“这两个分子能结合”而把机器学习引入靶点筛选、再用分子动力学模拟验证复合物稳定性整条链路才站得住。适合已经跑过 TCMSP、STRING、AutoDock 的人进阶也适合刚接触论文复现、想找一个端到端项目练手的人。下面按我实际复现的顺序把每一步的命令、参数和翻车点讲清楚。2. 网络药理学打底从 CQD 成分到血吸虫病靶点的交集网络药理学是整条链路的入口它的产出质量直接决定后面机器学习和分子对接有没有意义。这一步的核心逻辑是先拿到 CQD 的化学成分再预测这些成分的作用靶点同时收集血吸虫病的已知靶点最后取交集。听起来简单但成分筛选阈值和靶点来源选错后面全是白干。2.1 成分与靶点数据的获取和过滤常见做法是从 TCMSP、BATMAN-TCM、SwissTargetPrediction 三个库交叉取。TCMSP 的优势是自带 OB口服生物利用度和 DL类药性两个筛选指标业界默认阈值是 OB ≥ 30%、DL ≥ 0.18。但我要提醒一句这两个阈值不是铁律CQD 里某些含量低但活性强的成分会被误杀所以我会额外保留 OB 在 20%~30% 之间、DL ≥ 0.15 且已有文献报道有抗寄生虫活性的成分。import pandas as pd # 读取 TCMSP 导出的 CQD 成分表 herb pd.read_csv(cqd_ingredients.csv) # 标准筛选OB30 且 DL0.18 std herb[(herb[OB] 30) (herb[DL] 0.18)] # 放宽筛选OB 20-30 且 DL0.15作为补充候选 relaxed herb[(herb[OB] 20) (herb[OB] 30) (herb[DL] 0.15)] # 合并去重保留 MOL_ID 唯一 candidates pd.concat([std, relaxed]).drop_duplicates(subsetMOL_ID) print(f标准筛选 {len(std)} 个放宽后合计 {len(candidates)} 个) candidates.to_csv(cqd_candidates.csv, indexFalse)这段代码的关键在drop_duplicates(subsetMOL_ID)因为同一成分可能出现在多味药材里不去重会导致后面靶点频次统计虚高。参数上OB 和 DL 的列名要按你实际导出的表头改TCMSP 不同批次导出列名可能是ob、dl小写。跑完先看数量CQD 这类复方通常标准筛选后剩 80~150 个成分如果只剩二三十个说明阈值卡太死或者药材名没对齐。靶点预测这一步我一般用 SwissTargetPrediction 批量提交 SMILES取 Probability ≥ 0.1 的结果再用 UniProt 把靶点名统一成 Gene Symbol。血吸虫病靶点则从 GeneCards、OMIM、DisGeNET 三个库取并集关键词用 “schistosomiasis”“Schistosoma japonicum”“Schistosoma mansoni”。这里有个血泪经验GeneCards 的 Relevance Score 别一刀切我一般保留 score ≥ 1 的全部再人工核对一遍因为血吸虫病本身靶点数据就不多砍太狠交集会空。2.2 交集靶点与 PPI 网络的构建拿到两边靶点后取交集得到 CQD 抗血吸虫病的潜在靶点。接着把交集靶点丢进 STRING 做 PPI 网络物种选 Homo sapiens置信度设 0.4中等置信度导出 TSV。然后用 Cytoscape 或 Python 的 networkx 算度值degree度值排名前 15~20 的通常就是核心靶点。import networkx as nx import pandas as pd # STRING 导出的边表 edges pd.read_csv(string_interactions.tsv, sep\t) G nx.from_pandas_edgelist(edges, node1, node2) # 计算度值并排序 deg pd.Series(dict(G.degree())).sort_values(ascendingFalse) core_targets deg.head(20).index.tolist() print(核心靶点, core_targets) # 导出给 Cytoscape 做可视化 nx.write_gexf(G, ppi_network.gexf)nx.from_pandas_edgelist默认建无向图如果你要做有向分析得加create_usingnx.DiGraph()。度值排序后核心靶点往往是 AKT1、TNF、IL6、TP53 这类泛靶点这很正常但也是坑——泛靶点太多说明你的成分靶点预测太宽泛需要回头收紧 SwissTargetPrediction 的 Probability 阈值。这一步的产出交集靶点列表 核心靶点是下一章机器学习特征工程的输入所以文件命名和格式要规范我习惯存成intersection_targets.csv和core_targets.txt。3. 机器学习筛选关键靶点特征怎么造、模型怎么选网络药理学给出的交集靶点动辄上百个直接全丢去做分子对接既费算力又没重点。机器学习在这里的作用不是“预测疗效”而是给靶点做重要性排序把候选缩小到 10~20 个。很多人一上来就套随机森林但特征工程没做好模型给出的重要性排序就是玄学。3.1 特征矩阵的构造把靶点变成可学习的向量每个靶点要变成一行特征。我一般构造三类特征第一类是网络拓扑特征包括 degree、betweenness、closeness、clustering coefficient第二类是功能富集特征把靶点对应的 GO 和 KEGG 条目做 one-hot第三类是文献支持特征用 PubMed 检索 “靶点名 schistosomiasis” 的命中数取对数。标签怎么来这是复现里最容易卡住的地方。常见做法是用已知抗血吸虫药物如吡喹酮的靶点作为正样本从交集靶点里随机抽等量非药物靶点作负样本。import numpy as np import networkx as nx from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score # G 为上一章的 PPI 网络 feat [] for node in G.nodes(): feat.append({ target: node, degree: G.degree(node), betweenness: nx.betweenness_centrality(G)[node], closeness: nx.closeness_centrality(G)[node], clustering: nx.clustering(G)[node], }) feat_df pd.DataFrame(feat).set_index(target) # 假设 pos_targets 为吡喹酮已知靶点neg_targets 为随机负样本 feat_df[label] 0 feat_df.loc[feat_df.index.isin(pos_targets), label] 1 X feat_df.drop(columnslabel).values y feat_df[label].values rf RandomForestClassifier(n_estimators500, max_depth6, random_state42) scores cross_val_score(rf, X, y, cv5, scoringroc_auc) print(AUC:, scores.mean()) rf.fit(X, y) importance pd.Series(rf.feature_importances_, indexfeat_df.drop(columnslabel).columns) print(importance.sort_values(ascendingFalse))参数上n_estimators500是为了让重要性排序稳定max_depth6是防止过拟合因为样本量通常只有一两百。cross_val_score用 5 折AUC 能到 0.75 以上说明特征有区分度如果 AUC 在 0.5 附近别急着调模型先检查正负样本是不是有信息泄漏——比如负样本里混进了正样本的邻居节点。这里要区分一下机器学习做靶点排序和计算机视觉那套完全不同没有卷积、没有图像张量本质是表格数据的分类问题所以别被“机器学习”四个字吓到sklearn 足够。3.2 模型解释与关键靶点输出随机森林给的是特征重要性不是靶点重要性。要得到靶点排序我一般用两种方式一是把每个靶点的特征向量输入模型取预测概率作为打分二是用 SHAP 值解释每个靶点对预测的贡献。后者更稳但计算量大。import shap explainer shap.TreeExplainer(rf) shap_values explainer.shap_values(X) # 对每个靶点求 SHAP 绝对值之和作为重要性 target_score pd.Series(np.abs(shap_values[1]).sum(axis1), indexfeat_df.index) target_score target_score.sort_values(ascendingFalse) print(target_score.head(15)) target_score.to_csv(target_importance.csv)shap_values[1]取的是正类label1的贡献别取错索引。输出的 top 15 靶点就是后续分子对接的受体清单。这一步的坑在于如果正样本只有五六个SHAP 值波动会很大建议做 10 次不同随机种子的平均。另外靶点重要性高不代表它一定是真靶点只是统计上更相关最终还得靠分子对接和动力学验证。我一般会把 top 15 和网络药理学 degree top 20 取交集双重过滤后的靶点更可信。4. 分子对接把关键靶点和 CQD 成分对上分子对接要回答的是CQD 里的哪个成分能和机器学习筛出的哪个靶点结合得多紧。这一步的产出是结合能binding energy单位 kcal/mol数值越负结合越稳。业内经验阈值是 ≤ -5.0 kcal/mol 算有结合可能≤ -7.0 算结合较好。4.1 受体和配体的准备受体靶点蛋白从 PDB 下载优先选有共结晶配体的结构分辨率 ≤ 2.5 Å。下载后用 PyMOL 或 Discovery Studio 去水、去配体、加氢。配体CQD 成分从 PubChem 下载 SDF用 OpenBabel 转成 PDBQT。# 受体准备去水去配体加氢用 AutoDockTools 的 prepare_receptor prepare_receptor4.py -r receptor.pdb -o receptor.pdbqt -A hydrogens # 配体批量转换 for f in ligands/*.sdf; do obabel $f -O ${f%.sdf}.pdbqt --gen3d doneprepare_receptor4.py是 AutoDockTools 自带的脚本-A hydrogens表示加极性氢。配体转换时--gen3d会重新生成三维构象这一步很关键因为 PubChem 下载的 SDF 有时是二维的不转三维对接结果全是错的。批量转换后检查文件大小PDBQT 文件如果只有几百字节说明转换失败多半是 SDF 里没有正确的连接表。4.2 对接盒子与结合能计算对接盒子grid box要包住靶点的活性口袋。如果你不知道口袋位置用 PyMOL 的fpocket插件或 CASTp 预测。盒子中心设口袋几何中心尺寸一般 20×20×20 Å 起步大口袋可以放到 30。# 生成对接参数文件 cat dock.conf EOF receptor receptor.pdbqt ligand ligand.pdbqt center_x 12.5 center_y -3.2 center_z 8.7 size_x 22 size_y 22 size_z 22 exhaustiveness 32 num_modes 9 EOF # 用 AutoDock Vina 对接 vina --config dock.conf --out ligand_out.pdbqt --log ligand_log.txtexhaustiveness32是精度和耗时的平衡点默认 8 太快结果不稳调到 64 以上单次对接可能超过十分钟。num_modes9输出 9 个构象取结合能最低的那个。跑完看 log 里的 affinity如果所有成分结合能都在 -4 以上要么盒子没对准口袋要么受体加氢有问题。我一般会拿原配体做一次重对接redockingRMSD ≤ 2 Å 才说明参数可信这一步是很多人省掉的后悔药。5. 分子动力学模拟验证复合物到底稳不稳分子对接给的是静态快照分子动力学模拟MD才是看复合物在溶剂里跑一段时间后还稳不稳。这一步算力消耗最大但也是整条证据链最有说服力的一环。常见做法是用 GROMACS跑 100 ns看 RMSD、RMSF、回旋半径和氢键数量。5.1 拓扑生成与体系构建蛋白用pdb2gmx生成拓扑配体用 ACPYPE 或 CGenFF 生成力场参数。力场选 AMBER99SB-ILDN 配 GAFF2这是蛋白-小分子复合物的常用组合。# 蛋白拓扑 gmx pdb2gmx -f receptor.pdb -o protein.gro -water tip3p -ff amber99sb-ildn # 配体拓扑用 ACPYPE acpype -i ligand.pdb -b ligand -c bcc -n 0 # 合并蛋白和配体 gmx editconf -f protein.gro -o box.gro -c -d 1.0 -bt cubic gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top 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-d 1.0表示蛋白到盒子边缘至少 1.0 nm太近会导致周期性镜像相互作用。-bt cubic立方盒子最省事但如果是膜蛋白得换 triclinic。genion加离子中和体系别跳过带电体系不中和跑起来会报错。配体拓扑生成后要手动把ligand.itp的#include写进topol.top这一步漏了grompp会直接报 “no such moleculetype”。5.2 能量最小化、平衡与成品模拟MD 分三步能量最小化、NVT/NPT 平衡、成品模拟。每步的 mdp 文件参数不同别混用。# 能量最小化 gmx grompp -f minim.mdp -c ionized.gro -p topol.top -o em.tpr gmx mdrun -v -deffnm em # NVT 平衡 100 ps gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx mdrun -deffnm nvt # NPT 平衡 100 ps gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr gmx mdrun -deffnm npt # 成品模拟 100 ns gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr gmx mdrun -deffnm mdminim.mdp里nsteps50000emtol1000nvt.mdp和npt.mdp里nsteps50000100 ps步长 2 fstcoupl用 V-rescalepcoupl用 Parrinello-Rahman。成品模拟nsteps50000000对应 100 ns。跑完用gmx rms、gmx rmsf、gmx hbond分析。RMSD 在 20 ns 后稳定在 0.2~0.3 nm 说明复合物稳定如果一直往上飘要么对接构象不对要么力场参数有问题。这里有个坑配体力场参数如果没做 RESP 电荷拟合跑出来的构象可能完全散开所以 ACPYPE 的-c bcc别省。6. 复现避坑从数据到算力的五个真实翻车点6.1 成分靶点预测结果为空现象SwissTargetPrediction 返回的靶点列表为空或只有一两个。原因提交的 SMILES 格式不对或者成分分子量太大超出预测范围。解决先用 RDKit 检查 SMILES 合法性Chem.MolFromSmiles(smi)返回 None 就是格式错分子量超过 1000 的成分换用其他预测工具。6.2 机器学习 AUC 异常高现象交叉验证 AUC 达到 0.99。原因正负样本有信息泄漏比如负样本里混入了正样本的直接邻居。解决构造负样本时排除正样本的一阶邻居或者用更严格的划分方式按网络社区划分训练测试集。6.3 分子对接结合能全部为正现象所有成分的 affinity 都是正值。原因受体 PDBQT 没有正确加氢加电荷或者对接盒子中心偏离口袋。解决用prepare_receptor4.py重新处理受体检查-A hydrogens是否生效用 PyMOL 把盒子中心可视化确认在口袋内。6.4 MD 模拟报 “LINCS warnings”现象mdrun中途报 LINCS 约束警告甚至崩溃。原因时间步长太大超过 2 fs或体系里有原子重叠。解决能量最小化没收敛就进平衡回去把emtol调小到 100nsteps加到 100000时间步长保持 2 fs别贪心用 4 fs。6.5 分析结果和论文对不上现象自己跑的 RMSD 曲线和论文图差异大。原因论文可能用了不同的力场、不同的模拟时长或者初始构象不同。解决先确认论文的力场和时长如果没写清楚以自己体系的收敛性为准别硬凑。MD 本身有随机性不同随机种子结果会有差异跑三次取平均更稳。7. 把四步串成一条可复用的流水线复现完这一篇最有价值的不是某个具体结果而是把网络药理学、机器学习、分子对接、分子动力学模拟串成了一条可复用的流水线。我的习惯是把每一步的输入输出固定成文件契约cqd_candidates.csv→intersection_targets.csv→target_importance.csv→docking_results.csv→md_analysis.xvg。这样换一个复方、换一个疾病只要替换第一步的药材和疾病靶点后面全部能跑。进阶用法上我一般会加两个动作。一是把分子对接的 top 5 复合物都跑 MD而不是只跑一个因为对接打分最高的不一定 MD 最稳多跑几个对比才有说服力。二是用 MM-PBSA 算结合自由能比单纯的 RMSD 更能定量说明结合强度。下面这个命令是 GROMACS 里算 MM-PBSA 的常用流程# 用 gmx_MMPBSA 计算结合自由能 gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx -cg 1 13 -ct md.xtc -cp topol.top-cg 1 13指定蛋白和配体的组号组号要先用make_ndx确认。mmpbsa.in里startframe设 5000对应 10 ns 后endframe设 50000interval设 50这样取 900 帧算平均结果比单帧可靠得多。验证方法上我习惯做三件事一是重对接 RMSD 检查二是 MD 跑三次不同随机种子看 RMSD 是否收敛到同一水平三是把关键残基突变掉再对接看结合能是否显著下降。这三步做完结论基本能站住。说个我自己的教训。第一次复现这类论文时我图快网络药理学阈值卡得死机器学习正样本只用了三个分子对接盒子凭感觉设MD 只跑了 10 ns。结果写出来的东西自己都不信。后来老老实实把每一步的参数记在本子上哪个阈值改了什么后果哪个参数动了结果怎么变才慢慢摸到门道。这类多方法串联的复现快就是慢慢就是快。希望帮到你。本文还有配套的精品资源点击获取
