网络药理学 + 分子对接 + 分子动力学：论文验证部分怎么做

先看结论

让验证部分站得住的五件事

审稿人看验证部分时，主要核对下面五件事。前三件在计算开始前就要定下来。

WFCMS《网络药理学评价方法指南》的解读文章统计了 CNKI 收录的网络药理学论文（截至 2021 年 6 月）：对结果做过实验验证的不足 30%，其中约三分之一用了体内外实验，多数研究只通过文献推测结果。审稿人见过大量只做对接的稿件，对接已不被当作验证。

全流程

从化合物到论文图：每一步的输入、工具、输出和图

工具版本为 2026-10-10 的当前版本。第 1、9、10、12 步决定验证部分是否站得住，第 2 到 8 步审稿人主要检查阈值和版本。

数据库

常用数据库 2026 年现状与使用经验

下表按 2026-10-10 实际访问结果整理。方法部分对每个数据库写三项：版本或网址、检索日期、筛选条件。

WFCMS 指南解读文章检索了 CNKI 中以人参为对象的 13 篇网络药理学论文：即使用同一个数据库、同一组阈值，纳入的成分数也不一致，原因包括阈值改动和文献补充成分未说明。TCMSP 中 OB > 30% 且 DL > 0.18 的人参成分共 22 个，其中有预测靶点的 17 个。保存原始导出表并在附件中提供，是让别人能复现你的成分数和靶点数的最直接做法。

阈值

筛选阈值：常见设置、出处与写法

多数阈值没有统一标准。可靠的做法是选一个有出处的值，在看到结果前固定，并在方法中写明理由。

PPI 网络

靶点交集与核心靶点：度值中位数的做法和局限

常规做法是把成分靶点和疾病靶点统一为 HGNC 基因符号后取交集，交集基因提交 STRING，在 Cytoscape 中计算拓扑参数，取度值不低于中位数的节点作核心靶点。

这样得到的核心靶点高度同质。Diao 等 2026 年统计了 1,038 篇天然产物网络药理学研究（2023 年 10 月至 2024 年 6 月发表）：AKT1 出现在 48.8% 研究的核心靶点中，TNF 为 39.1%，EGFR 为 34.8%；只依赖数据库的研究中 AKT1 占 54.9%，整合了组学或实验数据的研究中为 19.6%。槲皮素在两类研究中分别占 71.9% 和 42.3%。

原因之一是 STRING 的文本挖掘证据偏向被研究得多的基因。下面的脚本先用 get_string_ids 核对每个基因符号映射到的蛋白，再从固定版本的 STRING API 取网络，标出只有文本挖掘证据的边，并按度值中位数列出核心靶点。本页实测（2026-10-10，STRING 12.5 API，Python 3 标准库）：输入 AKT1、TNF、IL6、PTGS2、ESR1、STAT3、MAPK1、CASP3、BCL2、VEGFA 10 个网药高频基因，required_score=700 时得到 33 条边，其中 5 条只有文本挖掘证据，AKT1–IL6 的 0.88 分全部来自文本挖掘；9 个有连接的节点度中位数为 8，有 5 个节点并列达到中位数。

同一次实测发现：STRING 12.0 和 12.5 中查询 VEGFA 返回的是 COL18A1（9606.ENSP00000352798），不报错；用 UniProt 号 P15692 或 Entrez ID 7422 查询无结果，STRING 11.5 中 VEGFA 仍能正常匹配。直接把基因列表提交到 network 接口时，VEGFA 会被当作 COL18A1 进入网络（required_score=400 时 COL18A1 带着 3 条边出现）。VEGFA 是网药论文中的高频核心靶点，提交前先核对映射结果，映射不一致的基因在方法中单独说明。

降低同质性的做法：用疾病转录组或蛋白组数据核对核心靶点在疾病中是否有变化；除度值外同时报告介数等其他指标；在附件中给出完整靶点列表，便于审稿人判断。

# 用法：python3 string_hubs.py overlap_genes.txt 700
# overlap_genes.txt：每行一个交集基因符号；第二个参数为 required_score（0-1000）
import sys, csv, io, statistics, urllib.parse, urllib.request
from collections import defaultdict

API = "https://version-12-5.string-db.org/api/tsv"


def post(method, **params):
    params.update({"species": 9606, "caller_identity": "np_validation_tutorial"})
    data = urllib.parse.urlencode(params).encode()
    with urllib.request.urlopen(f"{API}/{method}", data=data) as r:
        return list(csv.DictReader(io.StringIO(r.read().decode()), delimiter="\t"))


genes = [g.strip() for g in open(sys.argv[1]) if g.strip()]
score = int(sys.argv[2]) if len(sys.argv) > 2 else 400

# 1. 先核对基因符号映射：STRING 会把查不到的符号静默映射到别的蛋白
mapped = post("get_string_ids", identifiers="\r".join(genes), echo_query=1, limit=1)
hit = {m["queryItem"]: m for m in mapped}
for g in genes:
    if g not in hit:
        print(f"WARNING not found in STRING: {g}")
    elif hit[g]["preferredName"].upper() != g.upper():
        print(f"WARNING {g} mapped to {hit[g]['preferredName']} ({hit[g]['stringId']})")
ids = [hit[g]["stringId"] for g in genes if g in hit and hit[g]["preferredName"].upper() == g.upper()]

# 2. 只用映射一致的 STRING ID 取网络
rows = post("network", identifiers="\r".join(ids), required_score=score)

edges = {}
for x in rows:
    key = tuple(sorted((x["preferredName_A"], x["preferredName_B"])))
    other = max(float(x[c]) for c in ("nscore", "fscore", "pscore", "ascore", "escore", "dscore"))
    edges[key] = (float(x["score"]), float(x["tscore"]), other)

tm_only = [k for k, (s, t, o) in edges.items() if o == 0]
print(f"edges={len(edges)}  textmining_only={len(tm_only)}")
for a, b in tm_only:
    print(f"  textmining only: {a}-{b} score={edges[(a, b)][0]}")

deg = defaultdict(int)
for a, b in edges:
    deg[a] += 1
    deg[b] += 1
med = statistics.median(deg.values())
print(f"nodes={len(deg)}  median_degree={med}")
for g, d in sorted(deg.items(), key=lambda kv: -kv[1]):
    print(f"{g}\t{d}\t{'>=median' if d >= med else ''}")

python3 string_hubs.py overlap_genes.txt 700

富集分析

GO/KEGG 富集：写明背景基因集和多重检验校正

富集分析的两个常见漏洞是没有交代背景基因集和没有做多重检验校正。

Wijesooriya 等 2022 年审查了 186 篇使用富集分析的文章：95% 的过表达分析没有使用或没有说明合适的背景基因集，43% 没有做多重检验校正。在网药研究里，AGE-RAGE 和 PI3K-Akt 通路分别出现在 54.9% 和 54.7% 的研究中（Diao 等 2026）。

可行的做法是用全部成分靶点与疾病靶点的并集作背景，并在方法中写明背景大小、校正方法和 KEGG 访问日期。clusterProfiler 的 enrichGO() 中 ont 参数默认是 "MF"，不写只会得到分子功能；enrichKEGG() 默认在线下载当前 KEGG 数据，因此要记录访问日期。

本页实测（2026-10-10，R 4.5.3，clusterProfiler 4.18.4，org.Hs.eg.db 3.22.0，macOS arm64）：上面 10 个基因作交集基因，以它们在 STRING 12.5 中的互作伙伴共 191 个基因作背景，脚本运行通过。使用该背景时显著的 GO BP 条目为 170 个、KEGG 通路 20 个；不设 universe（默认以全部有注释的基因为背景）时分别为 1,493 个和 132 个；不写 ont 时得到的是 77 个 MF 条目。两种背景下排名第一的 KEGG 通路都是 AGE-RAGE signaling pathway in diabetic complications。示例背景只用于演示，实际分析用全部成分靶点与疾病靶点的并集。

# 依赖：BiocManager::install(c("clusterProfiler", "org.Hs.eg.db"))
library(clusterProfiler)
library(org.Hs.eg.db)

hits <- readLines("overlap_genes.txt")     # 交集基因（基因符号）
bg   <- readLines("background_genes.txt")  # 背景：全部成分靶点与疾病靶点的并集

hits_id <- bitr(hits, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID
bg_id   <- bitr(bg,   fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID

ego <- enrichGO(gene = hits_id, universe = bg_id, OrgDb = org.Hs.eg.db,
                keyType = "ENTREZID", ont = "BP",  # ont 默认是 "MF"
                pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2,
                readable = TRUE)
ekegg <- enrichKEGG(gene = hits_id, universe = bg_id, organism = "hsa",
                    keyType = "kegg", pAdjustMethod = "BH",
                    pvalueCutoff = 0.05, qvalueCutoff = 0.2)

write.csv(as.data.frame(ego),   "go_bp.csv", row.names = FALSE)
write.csv(as.data.frame(ekegg), "kegg.csv",  row.names = FALSE)
writeLines(c(paste("hits", length(hits_id)),
             paste("background", length(bg_id)),
             paste("KEGG accessed", Sys.Date())), "enrichment_meta.txt")

分子对接

分子对接验证：阈值的来历、阳性对照与重对接

分子对接在网药论文中的作用是为核心成分和核心靶点提供结构层面的合理性。对接打分不能证明结合。

Che 和 Zhang 2025 年审查了 35 篇 2024 年发表的网药论文：约 15 篇用 −5.0 kcal/mol 作亲和力阈值，其他从 −1.2 到 −7.0 不等，都没有可追溯的出处；不同对接软件的打分不能互相比较；35 篇中没有一篇设阳性对照；4 篇明确使用盲对接，另有 13 篇不同配体对接到蛋白不同区域，疑似盲对接。他们在 CASF-2016 的 285 个复合物上测试：Vina 指定位点对接的打分与实验亲和力相关系数 0.604，盲对接为 0.387；盲对接时 exhaustiveness 为 8 的中位 RMSD 为 3.37 Å，提高到 64 降到 2.21 Å，耗时约为 8 倍。

受体与配体准备、盒子设置和完整命令见 /learn/autodock-vina-docking。

1. 01

   选结构

   优先选人源、含共晶配体、结合口袋完整的晶体结构，写明 PDB ID 和分辨率。只有预测结构时，先检查口袋残基的 pLDDT（判读方法见 /learn/alphafold-results）。

2. 02

   定位点

   用共晶配体位置或文献报道的功能残基定义对接盒子。位点未知时先做口袋预测，再用第二个对接程序交叉核对。

3. 03

   重对接

   取出共晶配体后用同一套参数重新对接，最优姿态与晶体姿态的重原子 RMSD ≤ 2 Å，说明这套参数能复现已知结合模式（PoseBusters 采用同一标准）。

4. 04

   阳性对照

   在同一盒子、同一参数下对接该靶点的已知抑制剂或共晶配体。成分打分与阳性对照比较，接近或优于阳性对照的成分才进入 MD。

5. 05

   报告参数

   写明盒子中心和尺寸、exhaustiveness、随机种子和 Vina 版本，附 2D/3D 相互作用图和重对接叠合图。

分子动力学

分子动力学验证：副本、指标与 MM/PBSA

分子对接和分子动力学是先后两步：对接在近似刚性的受体上搜索配体姿态并打分，分子动力学在显式溶剂中检验这个姿态随时间能否保持。对接位点错了，MD 无法纠正。

时长没有统一标准。Communications Biology 2023 年发表的分子动力学可靠性清单要求的是收敛证据和每个条件至少 3 次独立模拟，没有规定固定纳秒数。判断平衡的经验做法是：先把每帧按蛋白骨架最小二乘叠合到参考结构，再算 RMSD，对曲线做滑动平均，从某个时间点起不再有整体上升或下降，才算结构平衡。

MM/PBSA 的已知局限（Genheden 和 Ryde 2015）：算出的能量通常偏大；两个配体亲和力相差小于 12 kJ/mol（约 2.9 kcal/mol）时排序常不可靠；在整个 PDBbind 上与实验的 r² 约 0.3，单个蛋白在 0 到 0.8 之间；熵项统计误差最大，大规模测试中加入熵项并不改善结果，因此常被省略；溶质介电常数取 2 到 4 时结果往往最好；单轨迹法通常比三轨迹法准确。论文中写明是否计算熵项和所用介电常数。

gmx_MMPBSA 当前版本 1.7.0。照抄旧教程要注意：&gb 中 extdiel 自 1.5.0 起默认 78.5（此前为 80.0）；用 Interaction Entropy 时要报告 σIE，σIE 大于约 3.6 kcal/mol 时结果不可靠。体系搭建、mdp 参数和命令见 /learn/gromacs-protein-ligand，RMSD、RMSF、氢键与 SASA 的计算和作图见 /learn/md-trajectory-analysis。

• 每个复合物至少 3 条独立模拟（不同初始速度随机种子），报告均值与标准差。
• 给出被分析性质已平衡的证据，说明平衡段和分析段如何划分。
• 至少一条从不同初始构象（例如对接第二名姿态）出发的模拟，说明结论不依赖初始构象。
• 提供体系设置表：盒子尺寸、总原子数、水分子数、盐浓度、力场与水模型、质子化状态。
• 写明软件版本，附初始坐标、输入参数文件和最终坐标。
• 阳性对照复合物用同样方案模拟。

审稿意见

审稿人常见质疑与应对做法

下表的质疑来自期刊审稿要求和方法学审查文章，应对做法可以在投稿前完成。

国内做法

TCMSP 靶点名转基因符号与 Cytoscape 作图：经验帖做法与实测

旧版 TCMSP 的靶点表给出的是蛋白名（target_name 列），做交集和富集前要转成基因符号。国内经验帖的通用做法是下载 UniProt 人源、reviewed 条目，按蛋白名匹配基因名（Excel VLOOKUP 或 R 脚本）。下面给出这一步的实测结果和需要人工处理的情况，以及 Cytoscape 导入时的常见问题。

# 用法：python3 tcmsp2symbol.py tcmsp_targets.txt > target_symbol.tsv
# tcmsp_targets.txt：每行一个 TCMSP 靶点名（target_name 列，去重后）
# 只依赖 Python 3 标准库；运行时从 UniProt 下载人源 reviewed 条目（约 2 万行）
import csv, io, re, sys, urllib.request

URL = ("https://rest.uniprot.org/uniprotkb/stream?query=organism_id:9606+AND+reviewed:true"
       "&fields=accession,protein_name,gene_primary&format=tsv")
rec, alt = {}, {}
with urllib.request.urlopen(URL, timeout=600) as r:
    for row in csv.DictReader(io.StringIO(r.read().decode()), delimiter="\t"):
        gene = row["Gene Names (primary)"]
        if not gene:
            continue
        pn = row["Protein names"].split(" [Cleaved")[0]   # 去掉切割产物说明
        rec.setdefault(re.sub(r"\s*\(.*$", "", pn).strip().lower(), set()).add(gene)  # 推荐名
        for name in re.findall(r"\(([^()]*)\)", pn):      # 括号中的别名
            alt.setdefault(name.strip().lower(), set()).add(gene)

names = [l.strip() for l in open(sys.argv[1], encoding="utf-8") if l.strip()]
n_rec = n_alt = 0
print("target_name\tsymbols\tmatch")
for n in dict.fromkeys(names):
    key = n.lower()
    if key in rec:
        genes, how = rec[key], "recommended"; n_rec += 1
    elif key in alt:
        genes, how = alt[key], "alternative"; n_alt += 1   # 别名匹配：逐条人工核对
    else:
        genes, how = set(), "none"                          # 未匹配：人工查 UniProt
    flag = "ambiguous" if len(genes) > 1 else how
    print(f"{n}\t{'/'.join(sorted(genes))}\t{flag}")
print(f"targets={len(dict.fromkeys(names))} recommended={n_rec} alternative={n_alt} "
      f"unmatched={len(dict.fromkeys(names)) - n_rec - n_alt}", file=sys.stderr)

交给 Agent

把这套流程交给 Agent

计算部分的步骤固定、文件多、跑得久，适合交给科学智能体执行，判断和实验仍由你负责。

一句话指令示例：「以黄芩为研究对象、溃疡性结肠炎为疾病：从 HERB 2.0 和 BATMAN-TCM 2.0 收集成分与靶点，按 OB ≥ 30%、DL ≥ 0.18 和 SwissTargetPrediction 前 15 位筛选；疾病靶点用 OMIM、DisGeNET curated 和 GEO 差异基因；构建 STRING 12.5 网络并标出只有文本挖掘证据的边；以全部成分靶点与疾病靶点的并集为背景做 GO/KEGG 富集；对前 5 个核心靶点做定位点对接、重对接和阳性对照；对得分最高的 3 个复合物和阳性对照各跑 3 个 100 ns 副本，做 MM/GBSA 并出图。」

1. 01

   检索与记录

   访问各数据库，保存原始导出表，记录版本和检索日期。

2. 02

   网络与富集

   运行 Python 和 R 脚本完成交集、PPI、核心靶点和富集分析，输出表格和图。

3. 03

   对接

   下载 PDB 结构，准备受体和配体，完成重对接、成分对接和阳性对照对接，输出打分表和相互作用图。

4. 04

   分子动力学

   在预装 GROMACS、AmberTools、ACPYPE 的环境中搭建体系，按需租用 GPU 跑多副本，完成轨迹分析和 gmx_MMPBSA 计算。

5. 05

   对抗审阅

   检查任务是否真正完成，例如重对接 RMSD 是否 ≤ 2 Å、每个复合物的副本是否齐全、图中是否标注均值和标准差。

• 工作区保留：成分表、靶点表、网络 tsv、富集结果、对接日志与姿态、MD 输入文件与分析结果、全部脚本和图。
• 你仍需核对：核心成分是否确实存在于所用药材及其含量、阈值是否符合你的研究设计、PDB 结构和阳性对照的选择是否合适。
• 实验验证的设计和湿实验由你完成。