方法文章

可重复的计算工作流程用于药物发现:标准化网络药理学与分子对接分析

462 次观看

DOI:

10.3791/70171

2026年4月24日

* These authors contributed equally

本文内容

摘要

本标准化方案将网络药理学与分子对接结合分子动力学(MD)模拟,用于药物发现。该方案建立了定量筛选标准和可重复的实验步骤,适用于利用公共数据集进行多靶点药物筛选,并提高结果的可靠性。

摘要

网络药理学与分子对接技术广泛应用于药物发现领域,但零散的工作流程和不一致的操作常导致研究结果的可重复性降低。本文介绍了一种将上述方法整合为可重复框架的标准化操作流程,用于药物筛选与机制探索,该流程分为三个连续阶段:数据准备、计算分析和验证。在准备阶段,通过吸收、分布、代谢、排泄和毒性(ADMET)标准(包括口服生物利用度、类药性及毒性预测)对公共数据库中的化合物库进行筛选;同时通过靶点预测并整合疾病相关数据库,全面识别药物-疾病相互作用的候选靶点。在计算分析阶段,对交集靶点进行基因本体(GO)和京都基因与基因组百科全书(KEGG)富集分析,并结合蛋白质-蛋白质相互作用网络的拓扑学分析以确定核心靶点;分子对接采用两种具有不同优势的标准化可选策略进行配置。两步渐进策略首先使用AutoDock Vina对化合物库进行高通量初筛,随后利用YASARA进行精确的再对接,该方法可消除高通量筛选中的假阳性结果,并生成与后续YASARA分子动力学(MD)模拟天然兼容的蛋白-配体复合物,避免因跨软件格式转换引起的结构偏差。单步策略则完全通过YASARA独立完成全部对接过程,操作流程更简洁,实验效率更高,适用于特定研究目标。在验证阶段,采用标准化的分子动力学模拟,通过均方根偏差(RMSD)和均方根波动(RMSF)等核心指标评估配体-蛋白复合物的稳定性。该统一且可重复的技术流程提升了网络药理学与对接研究的可靠性,有助于计算药物发现领域的跨研究比较。

引言

网络药理学是一种从整体网络视角解析药物与生物体之间相互作用模式的研究方法1。通过构建并分析涵盖“药物-成分-靶点-疾病-生物通路”的相互作用网络,定量识别药物发挥疗效的关键分子、核心通路及协同作用机制。该分析框架符合国际公认的网络药理学标准,强调多组学数据整合与拓扑网络分析,最终阐明药物的整体治疗效应,预测潜在不良反应,或为新药研发提供系统性指导。该策略已成功应用于多个领域,包括解析非甾体抗炎药(NSAIDs)抗新型冠状病毒(COVID-19)作用靶点相关信号通路,以及探索治疗骨肉瘤和2型糖尿病等疾病的作用机制2,3,4。这一优势推动了现代药理学从单一靶点向多靶点药物发现的范式转变。

分子对接是一种计算模拟技术,通过算法建模评估小分子化合物与生物大分子靶标之间的空间匹配性和相互作用强度,从而预测最佳的结合构象5,6。Morris 等人开发了 AutoDock4 和 AutoDockTools4,这些工具支持受体选择性柔性的自动对接,基于分子动力学和分子几何学,通过计算分子间能量差异来评估结合的稳定性和亲和力7。AutoDock 等经典工具采用半柔性配体-受体建模方法,已成为此类模拟中评估结合亲和力的金标准7

在药物发现中,网络药理学筛选技术常用于鉴定关键成分,随后通过分子对接验证这些成分与靶标蛋白的结合能力,从而建立药物与疾病之间的关联8。然而,当前的应用普遍存在实验设计不一致、操作流程不规范以及缺乏验证步骤等问题。一项针对中药相关网络药理学研究的系统性评估已记录了这一普遍现象,指出不同数据库间的数据异质性以及实验验证不足导致研究结果难以重复9。这些局限性不仅削弱了研究成果的可重复性和可靠性,也阻碍了计算预测向实验验证的转化。与缺乏实验验证的孤立网络药理学研究或参数任意的碎片化对接流程不同,本方案将两种方法结合,并采用标准化阈值和逐步操作流程。该框架消除了主观参数选择,确保不同操作者均可获得一致结果,实现流程的可重复性。

因此,建立标准化的操作流程,通过网络药理学筛选识别潜在药物靶点和通路,并结合分子对接验证药物-靶点结合活性,可为研究药物作用机制及候选药物筛选提供实验依据10。该方案特别适用于利用公开可用的组学、化学和蛋白质结构数据库,对药物成分和合成小分子进行多靶点药物筛选;但不适用于尚无解析晶体结构的靶点,或缺乏已知疾病相关相互作用网络的单靶点孤儿药物筛选。

访问受限。请登录或开始试用以查看此内容。

方案

本方案仅涉及对公开数据库的计算分析,不涉及人类受试者、脊椎动物或生物组织的使用。本节所述的所有摘要工作流程均在图1中予以图示说明。

替代药物组分工作流程图;涉及分子对接和基因靶点识别。
图 1:工作流程概要。 绿色矩形代表替代药物组分,红色矩形代表疾病,黄色椭圆包含所使用的网站和软件,橙色矩形包含获得的文件或数据以及关键步骤,紫色菱形代表最终所需结果。 请点击此处查看此图的放大版本。

1. 药物成分及靶点的获取

  1. 使用化学名称作为关键词,在 PubChem 数据库(https://pubchem.ncbi.nlm.nih.gov/)中进行检索,以获取相应的 SMILES(Simplified Molecular-Input Line-Entry System)字符串。
  2. 访问 ADMETlab 3.0 网站(https://admetlab3.scbdd.com/),在“Services”选项卡下选择 ADMET Evaluation 功能,输入 SMILES 字符串,然后点击 SUBMIT 按钮。
  3. 根据吸收(Absorption)、分布(Distribution)、代谢(Metabolism)、排泄(Excretion)、毒性(Toxicity)、药物化学性质(Medicinal Chemistry)和毒发色团规则(Toxicophore Rules)等指标筛选 ADMET 结果,仅保留各项指标均满足预设阈值标准的化合物(表 1)。
  4. 访问 ProTox 3.0 网站(https://tox.charite.de/protox3/index.php?site=home),输入经筛选后的化合物 SMILES 字符串,选择 TOX PREDICTION 模块,勾选所有希望预测的毒性类型(如器官毒性、致癌性等),并运行预测。
  5. 根据 ProTox 3.0 的预测结果,剔除毒性预测值超过预设安全阈值的化合物(表 2)。
  6. 将通过 ADMET 和 ProTox 3.0 筛选的化合物整理成结构化的药物组分数据库(例如 Excel 或 CSV 格式),包含化合物名称、SMILES 字符串和筛选状态等字段。
  7. 访问 SwissTargetPrediction 网站(https://swisstargetprediction.ch/),在物种下拉菜单中选择 Homo sapiens,输入药物组分数据库中各组分的 SMILES 字符串,点击 Predict targets 按钮,并收集所有概率评分大于 0 的预测靶点。
  8. 访问 SEA(Similarity Ensemble Approach)网站(https://sea.bkslab.org/),输入上述用于靶点预测的相同 SMILES 字符串,并筛选结果,仅保留“Target Key”字段中以_Human结尾且 p 值低于 0.05 的条目。
  9. 将 SwissTargetPrediction 和 SEA 获得的靶点列表合并为一个药物作用靶点库。通过 Uniprot(https://www.uniprot.org/)将靶点名称去重并标准化为官方基因符号(例如,遵循 HGNC 指南)。
    注:药物作用靶点库可保存为 CSV 文件以供后续使用。

表1:ADMETlab 3.0 药物安全性筛选的阈值标准。 该表格总结了关键理化性质、ADME参数、代谢相互作用、毒性终点、毒性通路及毒理结构规则的推荐截断值与分类范围。根据概率值(< 0.3、0.3 - 0.7、> 0.7)或定量范围,预测结果被分为三个风险等级(低、中、高),可用于新药早期发现阶段对化合物安全谱的系统性评估。请点击此处下载该表格。

表2:ProTox-3.0 在药物发现中用于毒性预测的阈值标准。 本表总结了ProTox 3.0预测的关键毒性终点,重点关注药物早期发现阶段药物安全性评估中的关键参数。每个终点均提供二元分类结果(“有活性”或“无活性”)及相应的概率评分(0–1),“有活性”表示具有潜在毒性风险。应优先关注器官毒性(肝毒性、心脏毒性)、毒性终点(致癌性、致突变性、免疫毒性)以及CYP代谢酶抑制,因为这些是导致临床研发失败的主要原因。若多个终点均显示“有活性”,提示化合物具有广泛的毒性潜力,应考虑降低其研发优先级。急性毒性通过预测的LD50值和GHS分类进行评估,其中1–3类(< 300 mg/kg)被视为高毒性。概率评分反映了各项预测的置信度。请点击此处下载该表格。

2. 疾病靶点的获取

注意:在筛选数据库时,应统一目标基因的命名规范,以避免因命名差异导致遗漏。

  1. 访问五个疾病相关数据库:OMIM(https://www.omim.org/)、Disgenet(https://disgenet.com/)、TTD(https://ttd.idrblab.cn/)、GeneCards(https://www.genecards.org/)和 PharmGkb(https://www.pharmgkb.org/)。应用以下针对各数据库的筛选标准:对于 GeneCards,筛选相关性评分 ≥ 1.0 的条目;对于 DisGeNET,选择与目标疾病相关的条目;对于 PharmGKB,通过选择 Gene 选项,将结果限制为与基因相关的条目;对于 TTD,保留“Disease”列与目标疾病匹配的条目。
  2. 针对每个数据库,使用目标疾病的正式名称(例如,阿尔茨海默病)作为搜索关键词,检索所有相关靶点。
  3. 将来自五个数据库的靶点列表整合至一个电子表格中。通过比较各列表中的基因符号,去除重复的靶点。
  4. 通过 Uniprot 将所有剩余靶点名称标准化为官方基因符号,以解决命名不一致的问题。将标准化且去重后的列表保存为疾病靶点库(CSV 或 Excel 格式)。
    注意:疾病靶点库可与药物作用靶点库(步骤 1.9)一并保存,供后续步骤 3 使用。

3. 常见药物-疾病靶点的获取

  1. 访问 Venny 2.1.0 网络工具(https://bioinfogp.cnb.csic.es/tools/venny/)。将药物作用靶点库(步骤 1.9)和疾病靶点库(步骤 2.4)分别导入 Venny 2.1.0 的两个输入框中,生成显示两组靶点重叠情况的维恩图。
  2. 从维恩图结果中提取交集靶点。将这些靶点标记为药物-疾病共有靶点(潜在相互作用位点),并保存为 CSV 文件。

4. 蛋白质-蛋白质相互作用(PPI)网络的构建与核心靶点分析

  1. 访问 STRING 数据库(https://cn.string-db.org/)。在下拉菜单中选择Homo sapiens作为物种。
  2. 将常见药物-疾病靶点(步骤 3.2)导入 STRING 的输入框。将“最小必需相互作用评分”参数设置为高置信度(0.700),然后点击搜索以生成蛋白质-蛋白质相互作用(PPI)数据。将 PPI 数据导出为 TSV(制表符分隔值)文件。
  3. 打开已预装 CytoNCA 插件的 Cytoscape 软件。通过文件 > 导入 > 从文件导入网络菜单,将 PPI 的 TSV 文件导入 Cytoscape。
  4. 点击应用程序 > CytoNCA > 打开启动 CytoNCA 插件。选择五个用于核心靶点筛选的参考指标:中介性(Betweenness)、接近性(Closeness)、度值(Degree)、特征向量中心性(Eigenvector)和局部平均连通性(LAC)。
    注意: 所使用的五个关键拓扑学指标分别为:中介性(介数中心性,衡量某靶点出现在网络中所有最短路径上的频率)、接近性(接近中心性,反映某靶点到网络中其他所有靶点的平均最短路径长度)、度值(局部连接度,量化某靶点与其他靶点之间的直接相互作用数量)、特征向量中心性(综合考虑靶点自身的连接性及其所连接靶点的重要性)以及 LAC(局部平均连通性,评估某靶点直接邻近节点之间的连接密度)。
  5. 点击工具 > 分析网络菜单启动网络分析,然后点击确定。将分析结果导出为 CSV 表格。
  6. 计算上述五个指标的中位数,并保留达到或超过中位数的靶点。重复步骤 4.5 多次,直至剩余 10 至 20 个靶点。
  7. 根据度值(Degree)指标对剩余靶点进行排序(从高到低),初步选取排名前 10 的靶点作为核心基因。将核心基因列表保存为 CSV 文件。
    1. 为减少假阳性结果,并确保仅结构适宜的靶点进入对接研究,需进一步评估其结构可行性和可成药性:检查 PDB 数据库中是否存在可用的高分辨率晶体结构(≤ 2.5 Å),或评估是否可构建可靠的同源建模结构;使用结合口袋预测工具确认是否存在合适的结合位点;并结合文献或功能数据库交叉验证其在疾病通路中的已知相关性。
    2. 对于缺乏可用结构、无可成药结合口袋或与疾病无关的靶点,应降低其在对接研究中的优先级。但本步骤获得的完整核心靶点列表仍可用于 GO 和 KEGG 富集分析,因为这些分析无需结构信息。
      注意:步骤 4.6 和 4.7 中的基因数量可根据需要调整。通常在步骤 4.6 后保留 10 至 20 个靶点,建议在步骤 4.7 中至少保留 10 个核心基因,以确保有足够的数据量进行可靠的 GO 和 KEGG 富集分析,并获得一致的可视化趋势。

5. GO 和 KEGG 富集分析与可视化

注意:本部分阐明基因在细胞组分、功能及细胞内通路水平上的作用。

  1. 访问 DAVID 网络工具(https://davidbioinformatics.nih.gov/home.jsp)。选择Gene List作为输入类型,并将核心基因导入输入框中。
  2. 将 Identifier 设置为 OFFICIAL_GENE_SYMBOL,并在 Select species 中选择Homo sapiens。然后点击Submit List以上传核心基因。
  3. 进行 GO 富集分析时,选择GOTERM_BP_DIRECTGOTERM_CC_DIRECTGOTERM_MF_DIRECT类别。
  4. 进行 KEGG 富集分析时,选择KEGG_PATHWAY类别。将 GO 和 KEGG 分析的显著性阈值均设置为 p < 0.05。
  5. 点击Functional Annotation Chart生成富集结果。将 GO 和 KEGG 结果导出为 CSV 文件。使用 R Studio 软件及 ggplot2 绘制前 10 个富集项/通路的柱状图或气泡图。
    注意:显示的条目/通路数量可根据需要进行调整。

6. 使用 Autodock Vina 进行分子对接

注意:步骤6和步骤7均为分子对接步骤。步骤6使用AutoDock Vina 1.1.2软件,而步骤7使用YASARA 10.3.16。使用YASARA有助于后续进行YASARA分子动力学模拟。如果需要AutoDock Vina的对接结果,则YASARA中的对接结果应与AutoDock Vina的结果保持一致。这可避免因软件切换导致的差异,并确保分子动力学模拟验证结果的可靠性。具体方法如下:使用LigPlot+(版本2.3)打开步骤6.31的“result.pdb”结果,生成二维相互作用图,识别与配体相互作用的关键残基;然后在YASARA的对接步骤7.18中选择这些关键残基,并设置对接盒子大小以覆盖结合口袋,从而最大程度地保证Vina与YASARA之间对接位点的一致性。随后,在步骤7.19选择最优对接构象时,应确保配体与受体之间的关键相互作用残基与AutoDock Vina结果中识别出的残基一致。该一致性要求侧重于保留关键的相互作用模式,而非原子位置的完全对应;由于力场参数化和侧链柔性的差异,周边残基构象的微小变化是可预期的。只要与活性位点关键残基的重要相互作用得以保留,即可认为对接结果在交叉验证意义上具有一致性。若无需进行AutoDock Vina对接(步骤6),可直接执行步骤7。

  1. 通过在 PubChem 数据库中搜索相应的 SMILES 字符串(步骤 1.1),获取名为 ligand.sdf 的药物化合物的 SDF(结构数据文件)。
  2. 使用 Chem3D 软件打开 SDF 文件。在“计算”选项下,选择 MM2,然后点击 Minimize Energy,以对化合物结构进行自由能最小化。
  3. 通过选择 File > Save As 将最小化后的结构保存为 ligand.mol2 文件。从 RCSB PDB 数据库(https://www.rcsb.org/;通过 PDB ID 或基因名称搜索)获取核心基因蛋白受体的 PDB(蛋白质数据库)格式文件,命名为 receptor.pdb。
    1. 优先选择分辨率 ≤ 2.5 Å 且具有已解析结合位点的结构(如可用)。选择结构时,应检查条目是否完整(例如,是否包含所有预期结构域,是否存在大的未解析环区)、是否存在可能影响配体结合的突变,以及是否包含功能重要的辅因子(如血红素、金属离子)或共结晶配体。
    2. 对于已知具有寡聚组装形式的靶标,需考虑单体或多个亚基形式是否适用于研究问题;若二聚体或更高阶相互作用相关,可下载生物组装体。所选结构将在后续步骤中进一步处理,因此初步检查有助于避免后续问题。
  4. 使用 PyMOL 软件打开 receptor.pdb。在命令行中输入 remove organic 并按 Enter,以从蛋白结构中移除小分子配体。
    注意:如果使用共结晶配体定义结合位点,应首先记录该配体的三维中心坐标,然后在 PyMOL 命令行中输入 remove organic 并按 Enter 删除共结晶小分子;否则,直接运行 remove organic 命令以移除共结晶小分子。
  5. 在命令行中输入 remove solvent 并按 Enter,以从蛋白结构中移除自由水分子;使用命令 select metal_cofactor, resn [目标辅因子残基名称] 来识别功能关键的金属离子或辅因子(如 HEM、Zn2⁺、Mg2⁺),并确认其保留在结构中。
  6. 在 PyMOL 中点击 File > Export Molecule > Save,将清洗后的受体导出为 receptor_clean.pdb。
  7. 在 UCSF Chimera 1.19 中打开 receptor_clean.pdb。点击 Tools > Sequence > Sequence 显示序列,检查结合位点附近是否存在缺失的环区(缺失区域以红色轮廓框标示)。若存在缺失环区,通过在序列窗口菜单中选择 Structure > Modeller (loops/refinement) 重建环区,选择 non-terminal missing structure,设置适当的模型数量(例如 5 个),然后执行计算。完成后,选择最合理的模型。
  8. 在 Chimera 中优化结构。使用 Rotamers 工具(Dunbrack 库)对选定残基进行侧链优化,通过 Columns 菜单添加 Clashes 和 H‑Bonds 以进行评估,并选择冲突最少(0–1 为佳)且具有有利氢键的构象。然后使用 Dock Prep(AMBER ff14SB)添加氢原子并分配电荷。最后使用 Minimize Structure 工具进行能量最小化,通过选择主链原子(sel @ca,c,n,o)并反向选择,启用 Fixed atoms 以固定主链。通过选择 File > Save PDB 将处理后的结构保存为 receptor_optimized.pdb。
    注意:对于有序性良好的残基,跳过侧链优化。Dock Prep 会自动处理质子化状态。最小化过程中应固定主链。
  9. 在 PyMOL 中重新打开 receptor_optimized.pdb,并定义标准结合位点。如果存在共结晶配体,使用其坐标为中心设置网格:记录配体中心,然后使用 remove organic 命令将其移除。若无共结晶配体,可根据文献中已知的关键残基(例如 select binding_site, resi XXX-XXX)定义结合位点,或通过口袋预测工具目视识别潜在结合口袋以验证视觉判断。记录所定义位点的三维中心坐标(x/y/z),用于网格框设置。
    注意:此处记录的坐标用于居中 AutoDock Vina 网格。对于基于残基的定义,应计算所选残基的几何中心;对于通过视觉或预测工具识别的口袋,则采用空腔中心。在定义结合位点时,必须考虑预期的对接策略是靶向正构(活性)位点还是变构位点。对于正构靶向,结合位点应基于共结晶配体或文献报道的保守活性位点残基定义。对于变构靶向,可使用口袋预测工具识别潜在的变构位点,特别是对于已知存在变构调控的靶标。在缺乏先验信息的情况下,可通过全局对接并聚类预测的结合热点来辅助识别潜在的变构位点。这种灵活性使本方案可适用于正构和变构药物发现项目。
  10. 在 PyMOL 中点击 File > Export Molecule > Save,将最终优化的结构导出为 receptor.pdb。
  11. 在 AutoDock Tools 4.2.6 中点击 File > Read Molecule 打开 receptor.pdb。定义柔性残基。点击 Edit > Flexible Residues > Select Residues,选择在配体结合时可能发生构象变化的结合位点残基(选择 ≤ 10 个残基)。
    注意:此步骤允许选定的侧链在对接过程中移动,以考虑诱导契合效应。
  12. 将含柔性残基的受体保存为 PDB 文件。点击 File > Save,在下拉菜单中选择 Write PDB。在“Available PDB records”窗口中勾选 ATOM 和 CONECT,点击 ADD,然后点击 OK。将文件保存为 receptor.pdb。
    注意:此 PDB 文件包含柔性残基信息,将用于生成 PDBQT 文件。
  13. 准备大分子用于对接。点击 Grid > Macromolecule > Choose,选择 receptor.pdb 文件,然后点击 Select Molecule。点击 File > Save As 将受体保存为 PDBQT 文件,并命名为 receptor.pdbqt。
    注意:AutoDock Tools 会分配电荷和原子类型,将受体保存为 AutoDock 原生的 PDBQT 格式,以备后续生成网格框和进行对接计算。
  14. 点击 Ligand 菜单,选择 Input,然后点击 Open。选择 ligand.mol2 并点击 OK。再次点击 Ligand 菜单,选择 Torsions,然后点击 Detect Torsions。AutoDock Tools 将自动识别配体结构中的可旋转键(例如烷基链中的单键,不包括肽键的酰胺键)。
  15. 在 Torsion Selection 窗口中,验证检测到的可旋转键(保留所有有效的可旋转键,排除芳香环键等刚性键)。点击 Set 确认扭转定义,然后点击 Close
    注意:保留有效的可旋转键可确保配体在对接过程中能采取不同构象(柔性配体),而受体保持刚性——这是 AutoDock Vina 中半柔性对接的核心。
  16. 再次点击 Ligand 菜单,选择 Output,然后点击 Save as PDBQT。将文件命名为 ligand.pdbqt,并保存在与 receptor.pdbqt 相同的目录中。
  17. 点击 Display 菜单,选择 Secondary Structure。点击 Display Only,然后选择 Lines 并点击 Undisplay,以简化蛋白视图。
  18. 点击 Grid 菜单,选择 Grid Box。调整 x、y、z(中心坐标)和 Spacing(Å) 值,将网格框置于蛋白的活性位点上。
    注意:若结合位点未知,可使用口袋预测工具(如 CASTp、DoGSite)识别潜在结合口袋。覆盖整个蛋白会显著增加假阳性率和计算成本,不推荐。
  19. 点击 File > Close saving current,然后点击 Grid > Output > Save GPF,将网格框设置保存为 Grid.gpf。
  20. 使用文本编辑器打开 Grid.gpf,记录文件中的 gridcentre(x、y、z 值)和 npts(size x、y、z 值)。
  21. 创建一个名为 Config.txt 的新文本文件,并输入以下内容:
    receptor = receptor.pdbqt
    ligand = ligand.pdbqt
    center_x = [Grid.gpf 中的 gridcentre x 值]
    center_y = [Grid.gpf 中的 gridcentre y 值]
    center_z = [Grid.gpf 中的 gridcentre z 值]
    size_x = [Grid.gpf 中的 npts x 值]
    size_y = [Grid.gpf 中的 npts y 值]
    size_z = [Grid.gpf 中的 npts z 值]
    energy_range = 5
    num_modes = 10
    将括号中的文本替换为 Grid.gpf(步骤 6.19)中的值。
    注意:energy_range 参数应设置为相对于最优结合模型的最大允许能量差,单位为 kcal/mol。例如,设置为 5 表示当能量差达到 5 kcal/mol 时,AutoDock Vina 将终止计算。此外,num_modes 指定生成的结合模型数量,通常设置为 10。
  22. 将 vina_split.exe 和 vina.exe 文件放置在与 receptor.pdbqt、ligand.pdbqt 和 Config.txt 相同的目录中。
  23. 打开 Windows 系统控制台,使用 cd 命令导航至该目录(例如,cd C:\DockingFiles)。
  24. 输入以下命令并按 Enter:vina.exe --config config.txt --log log.txt --out output.pdbqt
  25. 等待对接完成(耗时因系统而异)。将生成两个文件:log.txt(对接结果)和 output.pdbqt(最低能量配体结构)。为确保可重复性,需使用不同的随机种子进行三次独立对接运行。若前几个构象之间的 RMSD < 1.0 Å,则表明结果一致。
    注意: 作为经验参考,AutoDock Vina 的结合能(kcal/mol)可解释为:≤ -7(高亲和力,潜在活性构象),-7 至 -5(中等亲和力),≥ -5(低亲和力)。这些阈值依赖于具体系统,应通过实验数据验证。
    1. 为评估特定靶标的对接准确性及区分能力,建议采用两种互补的验证方法。使用晶体学配体进行重新对接验证,以评估该方案能否重现实验观察到的结合模式,RMSD < 2.0 Å 作为标准接受准则。
    2. 使用公共基准数据集(如 DUD-E)进行富集分析,以评估该方案区分真实活性化合物与性质匹配的诱饵化合物的能力;这包括计算 ROC 曲线(提供分类性能的全局度量)和富集因子(如 EF1%),以量化高排名部分中活性化合物的富集程度。这些验证步骤共同有助于建立适当的亲和力截断值,并确保针对目标类别实现可靠的筛选性能。
  26. 打开 PyMOL 软件。点击 File > Open 导入 output.pdbqt 和 receptor.pdbqt。点击 File > Save As 将组合结构保存为 result.pdb。
  27. 点击 File > New Session 清除 PyMOL 工作区,然后重新打开 result.pdb 以可视化配体-蛋白复合物。

7. 使用 YASARA 进行分子对接

注意:此步骤为后续分子动力学(MD)模拟提供精确的重新对接和预处理,同时也是对第6步高通量初筛结果的逐步验证。第6步使用AutoDock Vina这一高通量虚拟筛选的金标准工具,从化合物库中快速筛选出具有优良结合亲和力的候选分子。本步骤采用YASARA进行对接,因其对接模块与YASARA MD模拟平台完全兼容,可避免因文件格式转换和软件切换导致的结构偏差,为后续MD模拟提供标准化的初始复合物结构。对于第6步中AutoDock Vina筛选出的所有候选分子,本步骤获得的对接结果(包括在活性口袋中的结合构象及关键氨基酸相互作用)必须与AutoDock Vina的结果一致,且结合亲和力的相对排序趋势需保持相同,方可进入MD模拟。由于两种软件的计算算法不同,其绝对对接评分不可直接比较。该一致性要求可消除因软件差异导致的假阳性结果,确保候选分子结合特性的稳定性,保障后续MD模拟验证的可靠性与逻辑连贯性。

  1. 使用 OpenBabel 将 ligand.sdf 文件转换为 ligand.pdb 文件。
    注意:此处 OpenBabel 仅用于格式转换。配体的分子动力学参数化将在后续步骤中由 YASARA 自动完成。
  2. 打开 YASARA 软件。点击 文件 > 加载,选择 ligand.pdb 以导入配体。点击 编辑 > 清理 > 全部,以去除配体的结构缺陷。
    注意:此步骤执行基本的几何结构清理。随后,YASARA 将使用其内置的 AutoSMILES 技术自动为配体分配力场参数,该技术应用通用 AMBER 力场(GAFF)和 AM1-BCC 电荷,以确保与用于蛋白质的 AMBER14 力场兼容。此参数化对于对接和分子动力学(MD)模拟中的精确能量计算至关重要。
  3. 点击 选项 > 默认 pH,选择合适的 pH 值(例如,生理条件下的 7.4),然后点击 确定
  4. 点击 对接 > 力场,设置对接所用的力场,以确保与后续 MD 模拟的参数一致性。
    注意:在 YASARA 10.3.16 中,AMBER14 是推荐用于此药物发现工作流程的力场,因其对蛋白质具有全面的参数覆盖,并与标准 MD 模拟协议完全兼容。对于标准蛋白质残基,参数将从力场内置模板中自动分配;对于小分子配体,YASARA 将使用其内置的 AutoSMILES 技术自动进行参数化,分配 GAFF(通用 AMBER 力场)原子类型和 AM1-BCC 电荷。这确保了蛋白质与配体参数之间的兼容性,从而在对接和 MD 模拟中实现精确的能量计算。可根据实际使用的 YASARA 版本及系统的具体特性选择更合适的力场。
  5. 点击 模拟器 > 定义模拟单元 > 围绕所有原子,以设置工作边界。点击 模拟器 > 单元边界 > 周期性,以启用周期性边界条件。
  6. 点击 选项 > 选择实验 > 能量最小化,然后点击 运行,以最小化配体的能量。
  7. 点击 文件 > 另存为,将文件命名为 ligand.pdb,点击 确定 以覆盖原始的配体 PDB 文件。点击 文件 > 新建 以清空工作区,然后点击 文件 > 加载,选择 receptor.pdb 文件。
  8. 对蛋白质受体重复步骤 7.2 至 7.7,将处理后的文件保存为新的 receptor.pdb 文件。
  9. 点击 文件 > 新建,然后点击 文件 > 加载,同时选择 ligand.pdbreceptor.pdb。重复步骤 7.3 至 7.5,设置 pH 值、定义模拟单元,并为复合物启用周期性边界。
  10. 点击 处理器 > 设置 CPU,选择要使用的 CPU 核心数。点击 处理器 > 设置 GPU,选择 GPU 设备以加速计算。
  11. 点击 文件 > 另存为 > YASARA 场景,将文件命名为 sce\nesult.sce(若 sce 文件夹不存在,请先创建),然后点击 确定
  12. 点击 选项 > 宏&动画 > 设置目标,选择 sce\nesult.sce,点击 确定。点击 选项 > 宏&动画 > 运行宏,选择 dock_run.mcr 宏文件,点击 确定
  13. 点击 模拟器 > 定义模拟单元 > 围绕选定原子,并重复步骤 7.5,然后点击 继续 以启动对接。
  14. 等待对接完成。将生成以 yob 为后缀的文件;name.log 文件包含结合能和接触受体残基。
    注意:为确保分子动力学模拟验证的合理性,应在 YASARA 中选择与 AutoDock Vina 对接结果一致的对接结果。

8. 分子动力学模拟

  1. 点击 文件 > 新方法 以清空工作区。然后单击 文件 > 加载 > YASARA 对象 并选择 result.yob.
  2. 在场景内容面板(右侧)中,展开所有 Mol 条目。单击 编辑 > 分离 > 对象,选择全部 摩尔 序列面板中的内容,然后单击 好的.
  3. 点击 编辑 > 加入 > 对象,选择全部 摩尔 除第一个和最后一个条目(配体)外的所有内容,并点击 好的. 选择第一个 摩尔 进入并点击 好的 再次连接以重新结合蛋白。
  4. 继续对组件重新编号。选择 重新编号 编辑并点击 物体这会产生两个部分:第一部分是蛋白质-受体复合物,第二部分是小分子配体。
  5. 点击 编辑 > 转移,然后点击 对象 从下拉列表中选择选项。在序列面板中,首先通过单击其对应条目选择小分子配体内容,然后通过单击其条目选择蛋白质受体内容并点击 好的 以确认选择配对。
  6. 在下一个弹出窗口中,勾选以“转移过程中固定屏幕上原子”开头的选项,然后单击 好的.
  7. 重复步骤 7.2 至 7.5,然后点击 模拟器 > 温度 并选择 298K. 点击 文件 > 保存为 > YASARA Scene,将文件命名为 sce\nesultrun.sce,然后单击 好的.
  8. 点击 文件 > 新 以清空工作区。然后单击 选项 > 大体&视频 > 设定目标,选择 sce esultrun.sce,然后单击 好的.
  9. 确保在步骤 7.4 中选择的力场也用于分子动力学(MD)模拟;md_run.mcr 宏通常会继承当前的力场设置。点击 选项 > 大体&视频 > 播放宏,选择 md_run.mcr 宏文件,然后单击 好的 开始分子动力学模拟。
  10. 对蛋白-配体复合物进行三次独立的分子动力学模拟(3 × 100 ns),每次使用不同的初始速度,并对三条轨迹进行统计分析,以确保结果的可靠性。运行过程中将生成 sim 格式的文件。例如,若每 100 ps 保存一次轨迹,则一次 100 ns 的模拟将生成 1000 个以 sim 为后缀的文件。
  11. 完成步骤 8.10 后,单击 选项 > 大体&视频 > 设定目标,选择 sce esultrun.sce 文件,然后单击 好的.
  12. 点击 选项 > 大体&视频 > 播放宏,选择 md_analyze.mcr, md_analyzebindenergy.mcr,以及 md_analyzeres.mcr 并单击 好的.
  13. 完成全部三项分析后,将生成相应的数据文件:result_run_analysis.tab、result_run_bindenergy.tab 和 result_run_analysisres.tab。
  14. 首先分析 result_run_analysis.tab 文件,该文件提供了10个核心参数:Energy(体系总能量)、Bond(键能)、Angle(键角能)、Dihedral(二面角能)、Planarity(平面性能量)、Coulomb(静电能)、VdW(范德华能)、CA(蛋白质Cα原子的RMSD)、Backbone(蛋白质主链RMSD)和 HeavyAtoms(重原子RMSD)。
  15. 提取时间(ns)列及相应的参数列,以评估体系是否达到能量平衡。通过观察初始10–20 ns后势能在一个较小波动范围内趋于稳定,确认体系的稳定性。通过监测Cα原子、蛋白质主链及重原子的均方根偏差(RMSD),评估构象稳定性。当这些RMSD值达到平台期时,认为模拟体系在结构上已达到稳定。
  16. 对于典型大小的蛋白质-配体复合物,Cα原子和主链RMSD值稳定在2.5 Å以下,同时重原子RMSD低于3.5 Å,可作为构象稳定性的支持性指标。但关键判断标准应为RMSD轨迹中是否存在明确的平台期,此为主要且必需的判据,而非仅严格依赖上述数值。
    注意:这些阈值是经验性的,应结合特定蛋白质的大小和柔性进行解读。判断收敛性的决定性指标是持续的平台期,表明结构已围绕一个稳定的构象集合达到稳定状态。
  17. 接下来,分析 result_run_bindenergy.tab 文件,该文件提供了配体与靶标在模拟轨迹过程中的结合能。计算整个模拟期间的平均结合能。在 YASARA 的 MM-PBSA 实现中,正值越大表示结合越强。适度强且稳定的相互作用通常表现为平均结合能为正值且足够大(具体数值依赖于体系,但可通过已知结合物或实验数据进行校准),同时标准差相对于平均值较小(例如,变异系数 < 50 - 60%,反映出模拟过程中波动有限。
    注意:本步骤中报告的结合能是使用严格的 MM-PBSA 方法计算的,而默认的 YASARA 结合能宏则采用更快的近似方法(BoundaryFast)。默认近似方法适用于快速筛选或相对比较,而 MM-PBSA 方法更适用于获得更精确的绝对结合自由能。正如作者在 YASARA 宏头部明确指出的:更正的结合能值表示更强的结合,负值并不表示无结合。因此,用户应将正值解读为结合更强的指示,其数值大小取决于具体的蛋白-配体体系。
  18. 最后,分析文件 result_run_analysisres.tab,该文件提供每个残基的数据,包括残基编号(Residue ID)、均方根偏差(RMSD)、主链均方根偏差(Backbone RMSD)、重原子均方根偏差(HeavyAtoms RMSD)和均方根涨落(RMSF)。分析应聚焦于已识别的稳定生产阶段。首先,确定位于靶标活性位点内的残基(例如,距离配体5 Å范围内的残基)。然后,利用这些数据评估模拟过程中这些单个活性位点残基的构象稳定性。
    注意:在蛋白质-配体复合物中,活性位点残基的稳定性可参考以下经验性指标:在稳定阶段,RMSF 值低于 1.0 Å,RMSD 波动在 1 - 1.5 Å 范围内,通常表明局部构象保持良好。若残基的 RMSF 超过 2.0 Å,可能提示其具有较高灵活性;此类残基应映射至三维结构上,以判断其是否对应功能相关的柔性区域(如环区或表面区域),或提示结合口袋内存在潜在的不稳定性。这些数值标准并非绝对准则,主要判断依据应为无显著的构象漂移,且需结合整个系统的收敛性进行综合评估。
  19. 数据文件整理完毕后,将整理好的数据导入 Prism 以生成相应的图表。

访问受限。请登录或开始试用以查看此内容。

结果

在对氯雷他定治疗过敏性鼻炎(AR)进行网络药理学分析之后,选择氯雷他定与PTGS2之间的相互作用作为代表性案例研究,以逐步展示分子对接和分子动力学(MD)模拟方案的应用。本示例旨在演示工作流程的执行过程和数据分析方法,而非对该特定相互作用提供生物学验证。为与实验数据进行定量评估,建议用户将该方案应用于公共数据库中已知结合亲和力的、特征明确的体系。

靶点识别与重叠分析
从公共数据库中总共获得127个与氯雷他定相关的靶点。通过多个疾病数据库收集了与过敏性鼻炎(AR)相关的靶点,共得到2,620个靶点。如维恩图(图2)所示,氯雷他定与过敏性鼻炎之间共有57个重叠靶点。这57个共同靶点被选为潜在的治疗靶点,用于后续的网络分析和富集研究。

访问受限。请登录或开始试用以查看此内容。

讨论

意义与关键步骤
本方案结合了网络药理学、分子对接和分子动力学模拟,相较于单一方法或双组合流程具有显著优势,有助于解决当前药物发现中存在的关键效率低下和可靠性不足的问题。整个流程依赖于三个关键步骤,以确保其可靠性,每一步均针对计算药物筛选中的核心挑战。首先,通过多数据库整合(例如,使用 PubChem 获取小分子结构,5 个数据库获取疾病靶点,GO 和 KEGG 用于生物学机制与通路注释)并结合 ADMET 性质筛选,避免了依赖单一数据库带来的偏倚,确保仅有具备体内潜力的成分得以进入后续研究4。其次,采用多种拓扑指标(介数、接近度、度值、特征向量和 LAC)对成分-靶点-疾病网络进行拓扑分析,优先识别驱动疾病发生发展的核心基因,而非治疗相关性较低的边缘靶点3,11。第三,通过 AutoDockTools 对靶标蛋白进行严格制备并精确定义结合位点,以避免非特异性对接。靶标蛋白的制备包括去除共结晶配体和水分子,Morris 等人已证明...

访问受限。请登录或开始试用以查看此内容。

披露

所有作者声明均无利益冲突。

致谢

国家重点研发计划&中国国家重点研发计划(2024YFC3506300,2024YFC3506301),国家中医药管理局高水平重点学科—中医体质医学(编号:zyyzdxk-2023251),国家自然科学基金面上项目(82204948),中国教育部基础与交叉学科突破计划(JYB2025XDXM612),湖北省重大科技专项(2023BCA005),湖北时珍实验室首席科学家研究项目(HSL2024SX0002)

访问受限。请登录或开始试用以查看此内容。

材料

本文使用的材料清单
姓名公司目录编号评论
ADMETlab 3.0中国科学院上海药物研究所N/A用于预测配体药物代谢动力学和毒理学特征的ADMET(吸收、分布、代谢、排泄、毒性)性质在线平台(URL: https://admetlab3.scbdd.com/)
AutoDock Tools(AutoDock 4)斯克里普斯研究所AutoDock 4.2.6用于分子对接模拟的软件套件;包含用于对接的 AutoDock 4,以及用于准备蛋白质和配体输入文件(添加氢原子、分配电荷、设置可旋转键)、定义对接网格和分析对接结果的 AutoDockTools(ADT)。
AutoDock Vina斯克里普斯研究所AutoDock Vina 1.1.2开源分子对接软件;用于预测小分子配体与蛋白质受体之间的结合亲和力及结合构象
Chem3D珀金埃尔默信息学Chem3D 2024分子建模软件;用于构建、优化和可视化小分子配体的三维结构
CytoscapeCytoscape 联盟(系统生物学研究所)Cytoscape 3.10.3用于可视化和分析生物网络的开源软件;适用于构建和编辑基因/蛋白质相互作用网络
DAVID(功能注释、可视化与集成发现数据库)美国过敏与传染病研究所(NIAID)N/A用于功能注释与富集分析的在线工具;适用于对目标基因进行 GO(Gene Ontology)和 KEGG(Kyoto Encyclopedia of Genes and Genomes)通路富集分析(URL: https://david.ncifcrf.gov/
DisGeNET 数据库巴塞罗那超级计算中心(BSC)N/A基因-疾病关联数据库;用于识别与特定疾病相关的基因(URL: https://disgenet.com/)
GeneCards 数据库魏茨曼科学研究所N/A人类基因综合数据库;用于获取全面的基因信息(如表达、功能、疾病关联)(URL: https://www.genecards.org/)
LigPlus欧洲分子生物学实验室-欧洲生物信息学研究所(EMBL-EBI)LigPlus 2.3从三维坐标文件自动生成二维蛋白质-配体相互作用图的软件 示意图展示了氢键、疏水相互作用以及结合位点的残基 注册时使用学术邮箱即可获取 https://www.ebi.ac.uk/thornton-srv/software/LigPlus/ .
OMIM 数据库约翰斯·霍普金斯大学医学院(与美国国家生物技术信息中心合作)N/A在线人类孟德尔遗传数据库;用于检索遗传性疾病及其相关基因的信息(URL: https://www.omim.org/
OpenBabelOpenBabel 开发团队N/A开源化学工具箱;用于在不同软件平台之间转换分子文件格式(例如,从 .mol2 转换为 .pdb)
PharmGKB 数据库斯坦福大学N/A药物基因组学知识库;用于检索基因-药物相互作用和药物基因组学变异信息(URL: https://www.pharmgkb.org/
PrismGraphPad SoftwarePrism 9用于科学绘图、数据分析(例如绘制结合能分布曲线、分析误差棒)以及生成适合发表的高质量图表。
ProTox 3.0查里特é - 大学ä德国柏林tsmedizinN/A用于预测小分子毒理学终点的在线工具;可用于评估候选配体的潜在毒性(URL: https://tox.charite.de/protox3/index.php?site=home
PubChem 数据库美国国家生物技术信息中心(NCBI)N/A化学信息公共数据库;用于检索小分子配体的二维/三维结构及理化性质(URL: https://pubchem.ncbi.nlm.nih.gov/
PyMOL施ödinger, LLCPyMOL 2.6.1分子可视化软件;用于查看、编辑及生成蛋白质-配体复合物的高质量图像
R StudioPosit, PBCRstudio 2025.09.1+401用于R语言编程的集成开发环境(IDE);适用于生物数据的统计分析以及GO/KEGG图谱的绘制
RCSB PDB 数据库结构生物信息学研究协作实验室(RCSB)N/A蛋白质结构数据库;用于以 PDB 格式检索蛋白质受体的三维结构(URL: https://www.rcsb.org/
SEA(相似性集合分析方法)斯克里普斯研究所N/A基于化学相似性的靶点预测在线工具;用于补充SwissTargetPrediction以确认配体靶点(URL: https://sea.bkslab.org/)
STRINGSTRING 联盟(EBI、SIB 等)N/A已知和预测的蛋白质-蛋白质相互作用数据库;用于构建基因/蛋白质相互作用网络(URL: https://string-db.org/)
SwissTargetPrediction瑞士生物信息学研究所(SIB)N/A用于预测小分子潜在蛋白质靶点的在线服务器;可用于鉴定配体的候选受体(URL: http://swisstargetprediction.ch/)
TTD 数据库中山大学药物发现与开发研究所(IDRBL)N/A治疗靶点数据库;用于检索已验证及潜在药物靶点的信息(URL: https://db.idrblab.net/ttd/)
UCSF Chimera加州大学旧金山分校生物计算、可视化与信息学资源中心(Resource for Biocomputing, Visualization, and Informatics, RBVI)UCSF Chimera 1.19分子可视化与分析软件;用于蛋白质结构准备,包括缺失环区重建(通过 Modeller 接口)、侧链优化(Dunbrack 旋转异构体库)、质子化状态调整,以及采用 AMBER ff14SB 力场进行能量最小化。1.19 版本(2025 年 3 月发布)修复了 PDB 结构获取功能 . 可免费用于非商业用途,网址为 https://www.cgl.ucsf.edu/chimera/ .
UniProt数据库UniProt 联盟(EBI,SIB,PIR)N/A蛋白质序列与功能的综合数据库;用于检索蛋白质序列、结构及功能注释(URL: https://www.uniprot.org/
Venny 2.1.0国家生物技术中心í西班牙国家研究委员会生物分子与细胞生物学研究所(CNB-CSIC)N/A用于生成维恩图的在线工具;可用于可视化基因集之间的重叠情况(例如来自不同数据库的靶基因)(URL: https://bioinfogp.cnb.csic.es/tools/venny/)
YASARAYASARA BiosciencesYASARA 10.3.16分子建模与模拟软件;用于分子对接(步骤 3.7)及后续的分子动力学模拟,以验证对接结果

参考文献

  1. Hopkins, A. L. Network pharmacology: The next paradigm in drug discovery. Nat Chem Biol. 4 (11), 682-690 (2008).
  2. An, W., et al. Mechanisms of rhizoma coptidis against type 2 diabetes mellitus explored by network pharmacology combined with molecular docking and experimental validation. Sci Rep. 11 (1), 20849(2021).
  3. Hu, M., et al. Use of network pharmacology and molecular docking to explore the mechanism of action of curcuma in the treatment of osteosarcoma. Sci Rep. 13 (1), 9569(2023).
  4. Oh, K. K., Adnan, M., Cho, D. H. Network pharmacology approach to decipher signaling pathways associated with target proteins of NSAIDs against COVID-19. Sci Rep. 11 (1), 9606(2021).
  5. Kuntz, I. D., Blaney, J. M., Oatley, S. J., Langridge, R., Ferrin, T. E. A geometric approach to macromolecule-ligand interactions. J Mol Biol. 161 (2), 269-288 (1982).
  6. Sahu, M. K., Nayak, A. K., Hailemeskel, B., Eyupoglu, O. E. Exploring recent updates on molecular docking: Types, method, application, limitation & future prospects. Int J Pharma Res Allied Sci. 13 (2), 24-40 (2024).
  7. Morris, G. M., et al. Autodock4 and autodocktools4: Automated docking with selective receptor flexibility. J Comput Chem. 30 (16), 2785-2791 (2009).
  8. Li, C., et al. Characterization of the molecular mechanisms underlying lurasidone-induced acute manic episodes in bipolar depression: A network pharmacology and molecular docking approach. CNS Neurosci Ther. 31 (4), e70383(2025).
  9. Lee, W. Y., et al. Evaluating current status of network pharmacology for herbal medicine focusing on identifying mechanisms and therapeutic effects. J Adv Res. 76, 799-815 (2025).
  10. Shahzadi, Z., et al. Network pharmacology and molecular docking: Combined computational approaches to explore the antihypertensive potential of Fabaceae species. Bioresour Bioprocess. 11 (1), 53(2024).
  11. Che, X., Zhang, L. Blind docking methods have been inappropriately used in most network pharmacology analysis. Front Pharmacol. 16, 1566772(2025).
  12. Ren, M., Ma, J., Qu, M. Network pharmacology integrated with molecular docking and molecular dynamics simulations to explore the mechanism of shaoyao gancao tang in the treatment of asthma and irritable bowel syndrome. Medicine .(Baltimore). 103 (50), e40929(2024).
  13. Jorgensen, W. L. The many roles of computation in drug discovery. Science. 303 (5665), 1813-1818 (2004).
  14. Gao, L., et al. Molecular dynamics simulation-driven focused virtual screening and experimental validation of fisetin as an inhibitor of Helicobacter pylori htra protease. Mol Divers. 29 (6), 6243-6258 (2025).
  15. Schaefer, M. H., Serrano, L., Andrade-Navarro, M. A. Correcting for the study bias associated with protein-protein interaction measurements reveals differences between protein degree distributions from different cancer types. Front Genet. 6, 260(2015).
  16. Richter, S., Fetzer, I., Thullner, M., Centler, F., Dittrich, P. Towards rule-based metabolic databases: A requirement analysis based on KEGG. Int J Data Min Bioinform. 13 (3), 289-319 (2015).
  17. Gu, S., et al. Benchmarking ai-powered docking methods from the perspective of virtual screening. Nat Machine Intell. 7 (3), 509-520 (2025).
  18. Zhang, P., et al. Network pharmacology: Towards the artificial intelligence-based precision traditional chinese medicine. Brief Bioinform. 25 (1), 1-12 (2023).

访问受限。请登录或开始试用以查看此内容。

重印与许可

申请许可以重复使用本 JoVE 文章的文本或图表

申请许可

标签

ADMET KEGG

相关文章