企业官网建设流程全解析

简介：这份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") # 标准筛选：OB>=30 且 DL>=0.18 std = herb[(herb["OB"] >= 30) & (herb["DL"] >= 0.18)] # 放宽筛选：OB 20-30 且 DL>=0.15，作为补充候选 relaxed = herb[(herb["OB"] >= 20) & (herb["OB"] < 30) & (herb["DL"] >= 0.15)] # 合并去重，保留 MOL_ID 唯一 candidates = pd.concat([std, relaxed]).drop_duplicates(subset="MOL_ID") print(f"标准筛选 {len(std)} 个，放宽后合计 {len(candidates)} 个") candidates.to_csv("cqd_candidates.csv", index=False)

这段代码的关键在drop_duplicates(subset="MOL_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(ascending=False) core_targets = deg.head(20).index.tolist() print("核心靶点：", core_targets) # 导出给 Cytoscape 做可视化 nx.write_gexf(G, "ppi_network.gexf")

nx.from_pandas_edgelist默认建无向图，如果你要做有向分析得加create_using=nx.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(columns="label").values y = feat_df["label"].values rf = RandomForestClassifier(n_estimators=500, max_depth=6, random_state=42) scores = cross_val_score(rf, X, y, cv=5, scoring="roc_auc") print("AUC:", scores.mean()) rf.fit(X, y) importance = pd.Series(rf.feature_importances_, index=feat_df.drop(columns="label").columns) print(importance.sort_values(ascending=False))

参数上，n_estimators=500是为了让重要性排序稳定，max_depth=6是防止过拟合，因为样本量通常只有一两百。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(axis=1), index=feat_df.index) target_score = target_score.sort_values(ascending=False) print(target_score.head(15)) target_score.to_csv("target_importance.csv")

shap_values[1]取的是正类（label=1）的贡献，别取错索引。输出的 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 done

prepare_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.txt

exhaustiveness=32是精度和耗时的平衡点，默认 8 太快结果不稳，调到 64 以上单次对接可能超过十分钟。num_modes=9输出 9 个构象，取结合能最低的那个。跑完看 log 里的 affinity，如果所有成分结合能都在 -4 以上，要么盒子没对准口袋，要么受体加氢有问题。我一般会拿原配体做一次重对接（redocking），RMSD ≤ 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 md

minim.mdp里nsteps=50000，emtol=1000；nvt.mdp和npt.mdp里nsteps=50000（100 ps，步长 2 fs），tcoupl用 V-rescale，pcoupl用 Parrinello-Rahman。成品模拟nsteps=50000000对应 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调小到 100，nsteps加到 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设 50000，interval设 50，这样取 900 帧算平均，结果比单帧可靠得多。

验证方法上，我习惯做三件事：一是重对接 RMSD 检查，二是 MD 跑三次不同随机种子看 RMSD 是否收敛到同一水平，三是把关键残基突变掉再对接，看结合能是否显著下降。这三步做完，结论基本能站住。

说个我自己的教训。第一次复现这类论文时，我图快，网络药理学阈值卡得死，机器学习正样本只用了三个，分子对接盒子凭感觉设，MD 只跑了 10 ns。结果写出来的东西自己都不信。后来老老实实把每一步的参数记在本子上，哪个阈值改了什么后果，哪个参数动了结果怎么变，才慢慢摸到门道。这类多方法串联的复现，快就是慢，慢就是快。希望帮到你。

本文还有配套的精品资源，点击获取