简介：这份资源围绕除风清脾汤治疗血吸虫病的机制研究展开，面向具备生物信息学与机器学习基础的中医药现代化科研人员。内容整合网络药理学、机器学习、分子对接与分子动力学模拟，完整呈现从TCMSP、UniProt获取成分与靶点，构建草药-靶点网络，到Venn图筛选共同靶点、PPI分析、GO与KEGG富集，再到LASSO、随机森林、SVM-RFE筛选关键靶点，并验证汉黄芩素、山奈酚、木犀草素、槲皮素与靶点相互作用的分析链路，尤其突出山奈酚与TP53稳定结合这一发现。资源包为1个PDF文件，约808KB，内含详细可运行代码及逐段解释，可作为复现论文与迁移方法的参考模板。目前已有90人学习，适合希望掌握中药复方多成分-多靶点-多通路研究全流程的读者。
1. 从一份能跑通的代码包说起：网络药理学+机器学习复现到底值不值得下
血吸虫病这块，传统中药复方"除风清脾汤"（CQD）到底怎么起效，光靠"清热解毒"四个字说不清楚。这份资源做的事情，就是把网络药理学、机器学习、分子对接和分子动力学模拟串成一条完整链路，从TCMSP里捞成分、从疾病库里捞靶点、用LASSO+随机森林+SVM-RFE三件套筛核心基因，最后落到汉黄芩素、山奈酚、木犀草素、槲皮素和TP53的结合能上。整套代码是Python写的，NetworkX、PyVis、scikit-learn、MDAnalysis、seaborn全用上了，跑完能出Venn图、PPI网络、GO/KEGG气泡图、对接打分条形图和雷达图。适合谁？做中医药现代化、天然产物机制、或者想拿一个完整多组学分析模板练手的研究生和一线科研人员。不适合谁？指望复制粘贴就发SCI的——里面模拟数据占了大头，真数据得自己从数据库拉。但作为流程骨架和参数参考，这份东西能帮你省掉至少两周的踩坑时间。
2. 数据准备与网络构建：从TCMSP到草药-靶点互作图
2.1 为什么先做成分筛选而不是直接跑机器学习
网络药理学最怕的就是"垃圾进垃圾出"。CQD里化合物几百个，不先按OB（口服生物利用度）≥30%和DL（药物相似性）≥0.18筛一遍，后面PPI网络会炸成一团毛线球。代码里get_cqd_components()用的是模拟数据，但真实场景下你得去TCMSP逐个查。常见做法是：TCMSP搜"除风清脾汤"里每味药，把OB、DL、靶点列全导出来，再用UniProt的ID mapping把蛋白名转成基因symbol。这一步不做，后面Venn图交集会少得可怜。
import pandas as pd import numpy as np import networkx as nx from pyvis.network import Network def get_cqd_components(): # 真实场景：从TCMSP导出后按OB>=30, DL>=0.18过滤 cqd_components = { 'MOL_ID': ['MOL001', 'MOL002', 'MOL003', 'MOL004'], 'Molecule_Name': ['wogonin', 'kaempferol', 'luteolin', 'quercetin'], 'OB': [46.23, 41.88, 36.16, 46.43], 'DL': [0.24, 0.24, 0.25, 0.28], 'Target': ['PTGS2,ESR1,NOS2', 'PTGS2,ESR1,CASP3', 'PTGS2,ESR1,NFKB1', 'PTGS2,ESR1,CASP3,NFKB1'] } return pd.DataFrame(cqd_components) def build_herb_target_network(components_df): G = nx.Graph() for _, row in components_df.iterrows(): component = row['Molecule_Name'] targets = row['Target'].split(',') G.add_node(component, type='compound', size=15, color='blue') for target in targets: G.add_node(target, type='target', size=10, color='red') G.add_edge(component, target, weight=1) # PyVis输出交互式HTML，适合放补充材料 nt = Network(height='750px', width='100%', bgcolor='#222222', font_color='white') nt.from_nx(G) nt.show('herb_target_network.html') return G逻辑说明：build_herb_target_network把化合物和靶点分成两类节点，边代表相互作用。参数上，size控制节点大小，color区分类型，weight目前都是1，真实数据里可以用结合概率或文献支持度加权。注意PyVis生成的HTML在Jupyter里直接能看，但导出到PDF会丢交互，投稿时建议用NetworkX的matplotlib后端重绘静态图。
2.2 Venn图找交集：别小看这一步的坑
identify_common_targets()用venn库画CQD靶点和血吸虫病靶点的交集。真实场景下，疾病靶点从GeneCards、DisGeNET、OMIM三个库合并去重，通常能拿到几百到上千个。CQD这边筛完可能就几十个。交集往往只有个位数到十几个——如果交集少于5个，后面机器学习根本没法做，得回头放宽OB/DL阈值或者换数据库。
from venn import venn import matplotlib.pyplot as plt def identify_common_targets(components_df, disease_targets): all_compound_targets = set() for targets in components_df['Target']: all_compound_targets.update(targets.split(',')) plt.figure(figsize=(8, 6)) venn({'CQD Targets': all_compound_targets, 'Schistosomiasis Targets': set(disease_targets)}) plt.title('Common Targets between CQD and Schistosomiasis') plt.savefig('venn_diagram.png', dpi=300) plt.close() common_targets = all_compound_targets.intersection(set(disease_targets)) return list(common_targets)参数说明：dpi=300是投稿底线，别用默认的72。venn库对超过3个集合的支持一般，如果后面要加"健康对照"之类的第三组，建议换matplotlib-venn或者UpSetPlot。血泪经验：GeneCards的"Relevance score"别全要，通常卡中位数以上，不然假阳性靶点会把交集撑大，后面富集分析全是泛泛的通路。
3. PPI网络与机器学习筛选：三算法投票怎么定关键靶点
3.1 PPI网络拓扑：度中心性和中介中心性看什么
build_ppi_network()从STRING数据库拿互作关系，真实操作是去STRING网站输基因列表，下载string_interactions.tsv，里面combined_score列就是置信度。代码里模拟了7条边，真实数据通常几百条。度中心性（degree centrality）高的节点是"枢纽蛋白"，中介中心性（betweenness centrality）高的节点是"桥梁蛋白"。两者都高的，基本就是核心靶点。
def build_ppi_network(common_targets): # 真实场景：从STRING下载tsv，筛选combined_score > 0.4 ppi_data = [ ('PTGS2', 'ESR1', 0.9), ('PTGS2', 'CASP3', 0.8), ('ESR1', 'CASP3', 0.7), ('ESR1', 'NFKB1', 0.85), ('CASP3', 'NFKB1', 0.75), ('PTGS2', 'NFKB1', 0.8), ('NOS2', 'PTGS2', 0.6) ] G = nx.Graph() for source, target, score in ppi_data: if source in common_targets and target in common_targets: G.add_edge(source, target, weight=score) degree_centrality = nx.degree_centrality(G) betweenness_centrality = nx.betweenness_centrality(G) plt.figure(figsize=(10, 8)) pos = nx.spring_layout(G, seed=42) # seed固定，保证可重复 nx.draw(G, pos, with_labels=True, node_size=[v * 3000 for v in degree_centrality.values()], node_color=list(betweenness_centrality.values()), cmap=plt.cm.Blues) plt.title('PPI Network of Common Targets') plt.savefig('ppi_network.png', dpi=300) plt.close() return G, degree_centrality, betweenness_centrality关键参数：combined_score > 0.4是STRING的常用阈值，低于这个假阳性飙升。spring_layout的seed必须固定，不然每次跑出来的图布局不一样，审稿人会怀疑你数据有问题。节点大小用度中心性映射，颜色用中介中心性映射，一张图两个维度，比分开画两张更省版面。
3.2 LASSO+RF+SVM-RFE：三算法取交集才是稳的
单用LASSO会偏向选相关性强但可能冗余的基因，单用随机森林对噪声敏感，单用SVM-RFE计算量大且对参数敏感。三个一起跑，取排名都靠前的，才是相对稳的。代码里machine_learning_feature_selection()用模拟的100样本×7靶点矩阵，真实场景下样本量往往不够——这是网络药理学做机器学习的通病。常见做法是：用GEO数据库找血吸虫病相关的表达谱，或者用TCGA的泛癌数据做迁移，但要在文章里说清楚局限性。
from sklearn.ensemble import RandomForestClassifier from sklearn.svm import SVC from sklearn.linear_model import LogisticRegression from sklearn.feature_selection import RFE def machine_learning_feature_selection(common_targets): # 真实场景：X来自GEO表达谱，y是疾病/对照标签 np.random.seed(42) X = np.random.rand(100, len(common_targets)) y = np.random.randint(0, 2, 100) # LASSO：L1正则化，系数为0的直接淘汰 lasso = LogisticRegression(penalty='l1', solver='liblinear', C=0.1) lasso.fit(X, y) lasso_coef = pd.DataFrame({ 'Target': common_targets, 'Lasso_Coef': lasso.coef_[0] }).sort_values('Lasso_Coef', ascending=False) # 随机森林：看feature_importances_ rf = RandomForestClassifier(n_estimators=500, random_state=42) rf.fit(X, y) rf_importance = pd.DataFrame({ 'Target': common_targets, 'RF_Importance': rf.feature_importances_ }).sort_values('RF_Importance', ascending=False) # SVM-RFE：递归消除，ranking_为1的是最终保留 svc = SVC(kernel="linear") rfe = RFE(estimator=svc, n_features_to_select=3) rfe.fit(X, y) svm_rfe = pd.DataFrame({ 'Target': common_targets, 'SVM_RFE_Rank': rfe.ranking_ }).sort_values('SVM_RFE_Rank') result = pd.merge(lasso_coef, rf_importance, on='Target') result = pd.merge(result, svm_rfe, on='Target') key_targets = result.sort_values( by=['Lasso_Coef', 'RF_Importance', 'SVM_RFE_Rank'] ).head(3)['Target'].tolist() return key_targets, result参数说明：C=0.1是LASSO的正则化强度，越小惩罚越狠，选出的基因越少。n_estimators=500是随机森林的树数量，低于100不稳定，高于1000收益递减。n_features_to_select=3是SVM-RFE最终保留的特征数，通常根据样本量定，样本<50时选3-5个，样本>200可以选10个。注意：三个算法的排序标准不一样，LASSO看系数绝对值，RF看重要性，SVM-RFE看ranking，直接按列排序取head(3)是简化做法，更严谨的是取三者交集再按综合排名。
4. 分子对接与动力学模拟：从打分到稳定性验证
4.1 分子对接：AutoDock Vina才是正主
代码里molecular_docking()用np.random.uniform(-10, -5)模拟打分，真实场景必须用AutoDock Vina或LeDock。流程是：从PubChem下载化合物SDF，用OpenBabel转PDBQT；从PDB下载TP53晶体结构（比如1TUP），用PyMOL去水去配体，加氢后转PDBQT；写config.txt指定grid box中心坐标和大小，跑vina。结合能低于-7 kcal/mol算强结合，低于-9算很强。
def molecular_docking(key_components, key_targets): # 真实场景：调用AutoDock Vina，此处模拟结果 docking_results = [] for comp in key_components: for target in key_targets: score = np.random.uniform(-10, -5) docking_results.append({ 'Component': comp, 'Target': target, 'Docking_Score': score }) return pd.DataFrame(docking_results)真实操作命令示例（bash）：
# 准备受体和配体 obabel tp53.pdb -O tp53.pdbqt -xr obabel wogonin.sdf -O wogonin.pdbqt # 运行Vina vina --receptor tp53.pdbqt --ligand wogonin.pdbqt \ --center_x 10.5 --center_y 20.3 --center_z 15.8 \ --size_x 20 --size_y 20 --size_z 20 \ --exhaustiveness 32 --out wogonin_tp53.pdbqt参数说明：--exhaustiveness 32是搜索彻底度，默认8，调到32更准但更慢。--size_x/y/z是搜索盒子大小，一般20Å够用，太大浪费算力，太小可能漏掉结合位点。--center_x/y/z是盒子中心，用PyMOL看活性口袋坐标。注意：Vina的打分函数对金属酶和辅因子处理不好，TP53如果有锌离子，得在pdbqt里保留。
4.2 分子动力学模拟：MDAnalysis看RMSD和氢键
对接给的是静态快照，动力学模拟才能看稳定性。代码里analyze_tp53_docking()画了结合能和相互作用雷达图，但真实MD要用GROMACS或AMBER跑100ns以上，然后用MDAnalysis分析RMSD、RMSF、氢键数量。山奈酚结合最稳定的结论，通常来自RMSD波动小于2Å且氢键数量在模拟过程中保持稳定。
import MDAnalysis as mda from MDAnalysis.analysis import rms, hbonds def analyze_md_stability(topology, trajectory): u = mda.Universe(topology, trajectory) # RMSD分析 R = rms.RMSD(u, u, select='backbone') R.run() rmsd_df = pd.DataFrame(R.rmsd, columns=['Frame', 'Time', 'RMSD']) # 氢键分析 h = hbonds.HydrogenBondAnalysis(u, 'protein', 'resname LIG', distance=3.5, angle=120.0) h.run() return rmsd_df, h.hbonds参数说明：distance=3.5是氢键距离 cutoff（Å），angle=120.0是角度 cutoff（度），这是MD分析的标准值。select='backbone'只算骨架原子，避免侧链噪声。RMSD前10ns通常算平衡过程，从10ns后取平均和标准差。如果RMSD一直漂移不收敛，说明模拟时间不够或者体系没平衡好，得加长到200ns。
5. 避坑与常见问题：复现时最容易翻车的五个地方
5.1 现象：Venn图交集只有2-3个靶点，机器学习跑不起来
原因：TCMSP的OB/DL阈值卡太死，或者疾病数据库用了太严格的筛选条件。CQD里很多成分OB在30%边缘，DL在0.18边缘，稍微一动交集就变。
解决：先把OB降到20%、DL降到0.15试一次，看交集是否到10个以上。如果还是少，换用SwissTargetPrediction补靶点，或者把疾病数据库的Relevance score阈值从中位数降到25分位。记住：交集少于5个，后面所有分析都是空中楼阁。
5.2 现象：PPI网络图节点重叠成一团，审稿人说看不清
原因：spring_layout的k参数默认值不适合节点多的网络，或者没设seed导致每次布局不一样。
解决：节点超过50个时，k=0.5/sqrt(n)手动调，或者换kamada_kawai_layout。导出时用dpi=300以上，节点标签字号至少8。如果还是挤，只画度中心性前30的节点，其余用文字补充。
5.3 现象：LASSO跑出来所有系数都是0
原因：C参数太小，正则化太强，把所有特征都惩罚没了。或者X矩阵没标准化，量纲差异大。
解决：先对X做StandardScaler，然后C从1开始试，逐步降到0.01，看非零系数数量。如果始终为0，说明特征和标签真的没关系，得回头检查y标签是不是随机生成的。
5.4 现象：分子对接结合能全是-5到-6，没有强结合
原因：受体没去水、没加氢、或者grid box没对准活性口袋。
解决：PyMOL里remove solvent、h_add，然后用show cavity或者参考原配体位置定盒子中心。如果TP53有锌离子，obabel转pdbqt时加-xr保留金属。结合能普遍弱，试试换Vina的--scoring vinardo，对金属酶更友好。
5.5 现象：MD模拟RMSD一直上升不收敛
原因：体系没平衡好，或者模拟时间太短。常见于膜蛋白或大复合物。
解决：先跑NVT和NPT平衡各500ps，用gmx energy检查温度和压力是否稳定。RMSD前10ns算平衡，如果20ns还在涨，加到100ns。实在不收敛，检查力场选择——TP53这种核蛋白用CHARMM36m通常比AMBER99SB好。
6. 进阶技巧：把模拟数据换成真数据的三个关键操作
第一，TCMSP数据别手动抄。用requests直接调TCMSP的API（如果有）或者用selenium模拟浏览器导出，几百个化合物手动抄必出错。我一般写个脚本把OB、DL、靶点列全抓下来，存成CSV，后面所有分析都从这个CSV读。
第二，STRING的PPI数据下载后，用combined_score > 0.4过滤，然后导入Cytoscape做可视化。Cytoscape的yFiles Organic Layout比NetworkX的spring layout好看十倍，而且能导出矢量图。关键靶点用cytoHubba插件的MCC算法算一遍，和Python的度中心性对一下，两者都排前10的才写进文章。
第三，分子对接别只跑一个构象。Vina的--num_modes 9输出9个结合构象，选打分最好的那个做MD。但要注意，打分最好的不一定是最稳定的，我遇到过打分-9.5的构象跑MD 20ns就解离了，反而-8.2的稳定结合。所以MD验证这一步不能省，而且至少跑3个独立轨迹，每个100ns，看结果是否可重复。
# 批量处理对接结果，选最优构象 import glob import subprocess def run_vina_batch(receptor, ligands_dir, output_dir): results = [] for lig in glob.glob(f'{ligands_dir}/*.pdbqt'): name = lig.split('/')[-1].replace('.pdbqt', '') out = f'{output_dir}/{name}_out.pdbqt' subprocess.run([ 'vina', '--receptor', receptor, '--ligand', lig, '--center_x', '10.5', '--center_y', '20.3', '--center_z', '15.8', '--size_x', '20', '--size_y', '20', '--size_z', '20', '--exhaustiveness', '32', '--num_modes', '9', '--out', out ], check=True) # 解析输出，取第一条MODEL的打分 with open(out) as f: for line in f: if line.startswith('REMARK VINA RESULT'): score = float(line.split()[3]) results.append({'Ligand': name, 'Score': score}) break return pd.DataFrame(results).sort_values('Score')逻辑说明：--num_modes 9让Vina输出9个构象，REMARK VINA RESULT行第一个数字是结合能。批量跑完用pandas排序，取每个配体的最优打分。注意：subprocess.run的check=True会在Vina报错时抛异常，方便定位问题。如果配体多，建议用multiprocessing并行，但别超过CPU核数，不然Vina会抢资源。
从那以后我每次做网络药理学复现，都强制走一遍"交集靶点≥10 → PPI度中心性前20 → 三算法交集≥3 → 对接打分≤-7 → MD RMSD收敛"这条检查链，任何一环不满足就回头调参数，绝不硬着头皮往下跑。希望帮到你。
本文还有配套的精品资源，点击获取
