企业官网建设流程全解析

简介：这份资源面向具备生物信息学、网络药理学或中医药研究背景的科研人员，以及希望复现中药复方机制研究的学者，提供除风清脾汤治疗血吸虫病机制论文的完整复现代码与解释。内容围绕网络药理学、机器学习、分子对接与分子动力学模拟展开，涵盖TCMSP与UniProt靶点识别、草药-靶点网络构建、Venn共同靶点筛选、PPI网络分析、GO与KEGG富集，以及LASSO、随机森林、SVM-RFE关键靶点筛选和汉黄芩素、山奈酚等成分与TP53、TNF、IL6的对接验证。资源包为1个PDF文件，约808KB，内含可运行代码与逐段说明，便于对照论文流程理解多成分-多靶点-多通路机制。目前已有141人学习，适合作为类似中药复方研究的分析模板与代码参考。

1. 除风清脾汤治血吸虫病：一篇论文复现为什么值得你花两周

血吸虫病至今仍是部分流行区的公共卫生负担，吡喹酮是主力药，但耐药株报道和肝纤维化后遗症让“中药复方辅助治疗”成了不少课题组盯上的方向。除风清脾汤（CQD）是一张眼科古籍方，近年被拿来试治血吸虫病的肝肠病变，问题是：一张方子十几味药、上百个成分，到底谁在起作用、打在哪个靶点上？光靠动物实验说不清。这篇论文的思路是把网络药理学、机器学习、分子对接和分子动力学模拟串成一条证据链——先筛成分和靶点，再用机器学习从公共数据里找疾病关键基因，取交集后做对接，最后用分子动力学模拟验证结合稳不稳。适合做中药机制研究、天然产物筛选、或者想系统跑一遍“干实验+湿实验前预测”流程的研究生和一线科研人员。下面我按自己复现时的顺序，把每一步的命令、参数和翻车点讲清楚。

2. 网络药理学打底：从 TCMSP 到交集靶点的完整链路

网络药理学是这条链的起点，它的产出质量直接决定后面机器学习、分子对接有没有意义。很多人一上来就打开 TCMSP 搜药，结果成分表里一堆 OB、DL 不达标的化合物也照单全收，后面靶点爆炸，交集出来几百个基因，对接根本做不完。我一般会先把筛选阈值卡死，再谈别的。

2.1 成分筛选阈值怎么定：OB≥30% 与 DL≥0.18 的取舍

TCMSP 是中药系统药理学数据库，CQD 里的每味药（除风清脾汤常见组成如防风、荆芥、连翘、黄芩等，具体以论文原文为准）都要单独检索。核心两个参数：口服生物利用度 OB 和类药性 DL。行业默认门槛是 OB≥30%、DL≥0.18，这是 Lipinski 规则和大量中药研究沉淀下来的经验值。

但这里有个玄学：有些成分 OB 只有 25% 却是已知活性成分，一刀切会漏掉。我的做法是主表用 30/0.18，同时把 OB 在 20~30 之间、DL 达标的成分单独列一张“备选表”，后面如果交集靶点太少再回补。别小看这一步，血吸虫病本身靶点数据就偏少，交集太稀会让整条链断掉。

筛选完还要做名称标准化。TCMSP 给的靶点名是蛋白名（如 Prostaglandin G/H synthase 2），而 GeneCards、OMIM 给的是基因符号（PTGS2）。不统一，后面取交集全是空。

2.2 靶点预测与疾病基因取交集：一份能跑的 Python 脚本

成分靶点从 TCMSP 的 related targets 拿，疾病靶点从 GeneCards、OMIM、DisGeNET 搜“schistosomiasis”“Schistosoma japonicum”等关键词。两边都拿到后，用 UniProt 把蛋白名转成基因符号，再取交集。下面是我常用的处理脚本。

import pandas as pd # 读取成分靶点（TCMSP 导出，列名：molecule_name, target_name） comp = pd.read_csv("cqd_compound_targets.csv") # 读取疾病靶点（GeneCards 导出，列名：GeneSymbol） disease = pd.read_csv("schistosomiasis_genes.csv") # 蛋白名 -> 基因符号映射表，从 UniProt 下载后整理成两列 uniprot = pd.read_csv("uniprot_map.csv") # columns: protein_name, gene_symbol # 统一大小写、去空格，避免匹配失败 comp["target_name"] = comp["target_name"].str.strip().str.lower() uniprot["protein_name"] = uniprot["protein_name"].str.strip().str.lower() merged = comp.merge(uniprot, left_on="target_name", right_on="protein_name", how="left") # 映射不上的先留空，人工核查，不要直接丢 unmapped = merged[merged["gene_symbol"].isna()] print(f"未映射靶点数: {len(unmapped)}") comp_genes = set(merged["gene_symbol"].dropna().str.upper()) disease_genes = set(disease["GeneSymbol"].str.strip().str.upper()) intersection = comp_genes & disease_genes print(f"成分靶点 {len(comp_genes)} 个，疾病靶点 {len(disease_genes)} 个，交集 {len(intersection)} 个") pd.Series(sorted(intersection)).to_csv("intersection_targets.csv", index=False, header=["gene"])

逻辑说明：先做名称标准化再 merge，是因为 TCMSP 和 UniProt 的大小写、空格经常不一致，直接 merge 会大面积匹配失败。how="left"保留未映射的记录，方便你回头查是哪个靶点没转成功，而不是悄悄丢掉。参数上，intersection数量是关键指标——太少（比如低于 20）说明阈值太严或疾病关键词太窄，太多（超过 200）说明 DL 没卡住，两种情况都要回头调。

提示：GeneCards 免费版一次导出有数量限制，血吸虫病相关基因建议用多个关键词分别搜再合并去重，别只搜一个词。

2.3 PPI 网络与核心靶点：STRING 导出后怎么读

交集靶点丢进 STRING 做 PPI 网络，物种选 Homo sapiens，置信度 0.4 起步（论文里常见 0.4 或 0.7，0.7 更严但可能把网络打散）。导出 TSV 后用 Cytoscape 打开，跑 CytoNCA 或 cytoHubba 算 degree、betweenness。核心靶点一般取 degree 排名前 10~15 个，这些就是后面分子对接的受体候选。

这里容易翻车的是：STRING 默认会把你的基因名做一次自动校正，有些别名会被合并或丢弃。导出后一定核对节点数和你输入的交集数是否一致，少了就说明有基因没被识别，得手动换成官方 symbol。

3. 机器学习找疾病关键基因：别把分类器当万能钥匙

机器学习这一环的作用是从表达谱数据里找出真正区分“患病 vs 正常”的基因，给网络药理学交集做二次过滤。很多人以为随便上个随机森林就完事，其实数据预处理和验证方式才是决定结果可信度的东西。

3.1 数据从哪来：GEO 数据集下载与探针注释

血吸虫病相关的表达谱一般在 GEO 上找，关键词组合“Schistosoma”“liver fibrosis”“S. japonicum”。下载 Series Matrix 文件和平台注释文件（GPL）。原始数据是探针 ID，必须用 GPL 注释转成基因符号，一个基因对应多个探针时取平均表达值。

import pandas as pd import numpy as np expr = pd.read_csv("GSE_matrix.txt", sep="\t", index_col=0) gpl = pd.read_csv("GPL_annotation.txt", sep="\t", comment="#", low_memory=False) # 找到探针ID和基因符号两列，列名各平台不同，需按实际调整 probe2gene = gpl[["ID", "Gene symbol"]].dropna() probe2gene.columns = ["probe", "gene"] expr = expr.reset_index().rename(columns={"index": "probe"}) merged = expr.merge(probe2gene, on="probe", how="inner") # 多探针对应同一基因：取表达均值 expr_gene = merged.drop(columns=["probe"]).groupby("gene").mean() # log2 转换，判断是否需要：最大值远大于 20 通常说明未取对数 if expr_gene.values.max() > 50: expr_gene = np.log2(expr_gene + 1) expr_gene.to_csv("expr_gene_matrix.csv")

逻辑说明：groupby("gene").mean()处理多探针，是标准做法，但要注意有些平台一个探针跨多个基因，这种行在注释里是“GENE1 /// GENE2”，merge 前得先拆分或丢弃。log2 判断用最大值阈值，是因为不同数据集预处理程度不同，硬套公式会出错。参数上，how="inner"保证只保留能注释上的探针，注释率低于 60% 就要考虑换平台文件。

3.2 差异分析 + 特征选择：LASSO 与随机森林怎么配合

差异分析用 limma（R）或 Python 的 scipy 做 t 检验，阈值一般 |log2FC| > 1 且 adj.P < 0.05。得到差异基因后，用 LASSO 回归做特征选择，再用随机森林算重要性，两者取交集基因作为最终关键基因。LASSO 的 alpha 用交叉验证选，随机森林的 n_estimators 我一般设 500，max_depth 不限制先看重要性分布。

from sklearn.linear_model import LassoCV from sklearn.ensemble import RandomForestClassifier from sklearn.preprocessing import StandardScaler import pandas as pd X = pd.read_csv("expr_gene_matrix.csv", index_col=0).T y = pd.read_csv("sample_labels.csv", index_col=0)["label"] # 0 正常 1 患病 scaler = StandardScaler() X_scaled = scaler.fit_transform(X) # LASSO 交叉验证选 alpha lasso = LassoCV(cv=5, max_iter=10000, random_state=42).fit(X_scaled, y) lasso_genes = X.columns[lasso.coef_ != 0] # 随机森林重要性 rf = RandomForestClassifier(n_estimators=500, random_state=42).fit(X_scaled, y) importance = pd.Series(rf.feature_importances_, index=X.columns).sort_values(ascending=False) rf_genes = importance.head(50).index final_genes = set(lasso_genes) & set(rf_genes) print(f"LASSO 选出 {len(lasso_genes)}，RF 前50 {len(rf_genes)}，交集 {len(final_genes)}")

逻辑说明：LASSO 擅长剔除共线性冗余，随机森林能捕捉非线性关系，两者交集比单用一种稳。random_state固定是为了结果可复现，这点在论文复现里是硬要求。参数上，cv=5是样本量不大时的常规选择，样本少于 30 建议改 LOOCV；n_estimators=500是精度和耗时的平衡点，再往上收益很小。

注意：如果样本量只有个位数，机器学习结果基本不可信，这时候应该老实说明局限，而不是硬凑一个“关键基因”列表。

3.3 把机器学习结果和网络药理学交集对齐

机器学习得到的关键基因，和第二章的网络药理学交集靶点再取一次交集，这才是最终进入分子对接的靶点。这一步是整个流程的收敛点，通常从几百个基因收敛到 10~30 个。如果交集为空，说明两个数据源方向不一致，要么疾病关键词搜错了，要么表达谱数据集选的不是血吸虫病模型，得回头查。

4. 分子对接：从 PDB 取受体到结合能判读

分子对接是把前面筛出的成分和靶点两两配对，预测结合模式和亲和力。这一步的坑集中在受体准备和参数设置上，很多人跑完一堆负分就以为成功了，其实构象根本没对上。

4.1 受体和配体准备：PDB 下载与去水去配体

受体结构从 PDB 下载，优先选有共结晶配体的、分辨率高的（< 2.5 Å）。下载后用 PyMOL 或 AutoDockTools 去水、去原配体、加氢。配体从 PubChem 下载 SDF，用 Open Babel 转成 PDBQT。

# 受体准备：去水去配体加氢（PyMOL 命令行） pymol -cq receptor.pdb -d "remove solvent; remove resn HOH; remove organic; h_add; save receptor_clean.pdb" # 配体批量转换 for f in ligands/*.sdf; do obabel "$f" -O "${f%.sdf}.pdbqt" --gen3d done

逻辑说明：remove organic去掉原配体，是因为原配体占着活性口袋，不去掉对接结果会假阳性偏高。--gen3d让 Open Babel 生成三维构象，二维 SDF 直接对接是无效的。参数上，加氢的 pH 默认 7.4，如果你的靶点活性中心有特殊质子化状态，得手动调。

4.2 AutoDock Vina 参数：exhaustiveness 与 box 怎么设

对接用 AutoDock Vina，核心参数是--exhaustiveness（搜索彻底度）和--center_x/y/z、--size_x/y/z（对接盒子）。盒子要包住活性口袋，中心点取原配体几何中心最稳。

vina --receptor receptor.pdbqt --ligand ligand.pdbqt \ --center_x 12.5 --center_y 8.3 --center_z 20.1 \ --size_x 22 --size_y 22 --size_z 22 \ --exhaustiveness 32 --num_modes 9 --out out.pdbqt

逻辑说明：exhaustiveness默认 8，做论文复现建议 32 起步，值越高搜索越充分但耗时线性增长。盒子尺寸 22 Å 是常见经验值，太小装不下配体会导致结合能虚高，太大则搜索空间爆炸、结果发散。num_modes 9输出 9 个构象，取打分最高的那个，但一定要看它是否落在口袋内。

4.3 结合能判读：-7 kcal/mol 不是金标准

结合能（affinity）越负越好，一般 < -5.0 kcal/mol 算有结合可能，< -7.0 算较强。但这是经验阈值，不是金标准。我见过 -9 的构象其实挂在蛋白表面，也见过 -6 的正好卡在催化残基上。所以判读要结合三件事：结合能、构象位置、关键残基相互作用（氢键、π-π 堆积）。用 PyMOL 或 Discovery Studio 可视化，看配体是否和已知活性残基形成氢键，这比单纯看分数靠谱得多。

5. 分子动力学模拟：GROMACS 跑通 100 ns 的关键设置

分子对接给的是静态快照，分子动力学模拟（MD）才是验证结合稳不稳的手段。GROMACS 是主流选择，流程是：拓扑生成 → 能量最小化 → NVT/NPT 平衡 → 成品模拟 → 分析 RMSD/RMSF/结合自由能。

5.1 拓扑与力场：pdb2gmx 选什么力场

蛋白用pdb2gmx生成拓扑，力场选 AMBER99SB-ILDN 或 CHARMM36，前者对小分子兼容性好。配体拓扑用 ACPYPE 或 CGenFF 生成，力场要和蛋白一致。

# 蛋白拓扑 gmx pdb2gmx -f receptor_clean.pdb -o receptor.gro -water tip3p -ff amber99sb-ildn # 配体拓扑（ACPYPE） acpype -i ligand.pdb -b ligand -c bcc -n 0

逻辑说明：-water tip3p是常规水模型，和 AMBER 力场配套。-c bcc让 ACPYPE 用 AM1-BCC 算电荷，比默认的 Gasteiger 准，但慢。参数上，配体电荷方法必须和蛋白力场兼容，混用会导致体系能量异常。

5.2 平衡三步走：NVT、NPT 和位置限制

能量最小化后，先 NVT（恒温）100 ps，再 NPT（恒压）100 ps，期间对蛋白重原子加位置限制，让溶剂和离子先弛豫。之后去掉限制跑成品模拟。

# NVT 平衡 gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx mdrun -deffnm nvt # NPT 平衡 gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -p topol.top -o npt.tpr gmx mdrun -deffnm npt

逻辑说明：-r指定位置限制参考结构，NVT/NPT 阶段必须加，否则蛋白在溶剂没平衡好时就会变形。nvt.mdp里define = -DPOSRES开启限制，tcoupl用 V-rescale，pcoupl用 Parrinello-Rahman。参数上，温度 300 K、压力 1 bar 是生理条件默认值。

5.3 100 ns 成品模拟与 RMSD 判稳

成品模拟一般跑 100 ns，dt = 0.002 ps，每 10 ps 存一帧。跑完用gmx rms算 RMSD，看是否在 20 ns 后收敛到平台。RMSD 波动小于 0.2 nm 且无持续漂移，说明体系稳定。

gmx rms -s md.tpr -f md.xtc -o rmsd.xvg -tu ns gmx rmsf -s md.tpr -f md.xtc -o rmsf.xvg -res

逻辑说明：-tu ns让横轴单位是纳秒，方便读图。RMSF 看每个残基的柔性，活性口袋残基 RMSF 低说明结合稳定。参数上，如果 RMSD 一直漂移不收敛，先查是不是盒子太小导致蛋白和自身镜像作用，或者力场/电荷有问题。

6. 复现避坑：五个让我重跑过的坑

6.1 靶点名称不统一，交集直接为空

现象：网络药理学交集跑出来 0 个基因。原因：TCMSP 用蛋白名，GeneCards 用基因符号，没做映射。解决：先过 UniProt 映射表，映射不上的手动核查，别用模糊匹配硬凑。

6.2 机器学习样本太少，结果不可复现

现象：换 random_state 后关键基因列表大变。原因：样本量小（< 20），模型过拟合。解决：用 LOOCV，报告稳定性，或者干脆放弃机器学习这一环，只用网络药理学+对接，别硬凑。

6.3 对接盒子设错，结合能全是假阳性

现象：所有配体结合能都 < -8，但构象都在蛋白表面。原因：盒子中心没对准活性口袋。解决：用原配体几何中心定 box，对接后必须可视化确认构象在口袋内。

6.4 力场不兼容，MD 体系炸开

现象：能量最小化后 MD 第一步就报错或蛋白飞散。原因：配体电荷方法和蛋白力场不匹配。解决：统一用 AM1-BCC 电荷 + AMBER 力场，或统一用 CGenFF + CHARMM36，别混。

6.5 RMSD 不收敛，误判为结合不稳

现象：100 ns RMSD 一直上升。原因：盒子太小、周期性镜像作用，或初始结构有冲突。解决：检查盒子边缘到蛋白最小距离 > 1.0 nm，能量最小化步数加到 50000 步。

7. 进阶：用 MM-PBSA 算结合自由能，把结论钉死

跑完 100 ns 只能说明“体系稳”，但配体和靶点到底结合多强，得用 MM-PBSA（或 MM-GBSA）算结合自由能。这是把分子动力学从“看起来稳”推进到“能量上有利”的关键一步，也是审稿人最爱追问的地方。

流程是：从成品模拟轨迹里抽 100~200 帧，用 gmx_MMPBSA 或 AMBER 的 MMPBSA.py 计算。输入需要拓扑、轨迹、以及一个索引文件把受体和配体分开。

# 用 gmx_MMPBSA 计算，需先准备 mmpbsa.in gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx \ -cg 1 13 -ct md.xtc -cp topol.top

mmpbsa.in里关键设置：startframe和endframe取轨迹后 80% 的帧（前 20 ns 还没平衡，算进去会污染结果），interval取 5 或 10 抽帧。-cg 1 13是受体和配体的组号，得先用gmx make_ndx建好索引。

结果看 ΔG_bind，负值越大结合越强。但要注意：MM-PBSA 的绝对值误差较大，一般只用来横向比较同一靶点不同配体的相对强弱，别拿它和实验 IC50 直接对标。我一般会同时算每个残基的能量分解（per-residue decomposition），找出贡献最大的几个残基，和对接里的关键残基对照——如果两套方法指向同一批残基，这个结论就基本钉死了。

还有个细节：熵项（-TΔS）算起来很贵且噪声大，常规做法是省略，只报焓贡献，然后在文里说明。如果审稿人要求，再用正态模式分析补，但那是另一个量级的工作量。

我自己复现这类论文最大的教训是：别急着跑对接和 MD，前面网络药理学的靶点映射和机器学习的数据清洗才是决定成败的地方，这两步偷懒，后面全是白跑。另外，任何一步的结果都要可视化确认，分数只是参考，眼睛看到的构象和残基相互作用才是硬证据。希望帮到你。

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