基于WGCNA和网络药理学探讨小檗碱治疗结直肠癌的作用机制

WGCNA and Network Pharmacology Reveal Targets and Mechanisms of Action of Berberine in Colorectal Cancer Treatment

摘要: 目的：探讨小檗碱治疗结直肠癌的潜在治疗靶点和分子机制。方法：使用Swiss Target Prediction、Comparative Toxicogenomics、Target Net和Binding数据库筛选出小檗碱的活性成分及基因靶点。同时，基于GEO数据库的GSE90627数据集进行WGCNA分析，以识别结直肠癌的相关靶点。通过绘制维恩图以获得疾病–药物基因集。利用STRING数据库构建了蛋白互作网络(PPI)。利用Cytoscape软件进行拓扑分析，得到关键聚类。对关键簇基因进行了功能富集分析。最后，通过在线分子对接网站进行虚拟验证。结果：经过一系列的生物信息学分析，获得了包含35个基因的关键簇。富集分析结果表明，关键簇基因主要富集在p53和FOXO信号通路中。使用PPI和WGCNA分析鉴定出4个最重要的生物学靶标：WD重复蛋白74 (WDR 74)、细胞周期蛋白依赖性激酶4 (CDK 4)、极光激酶A (AURKA)和细胞周期蛋白1 (CCNE 1)。分子对接结果表明，上述4个靶标均能与小檗碱稳定结合。结论：本研究利用生物信息学和网络药理学方法筛选和鉴定了小檗碱治疗结直肠癌的生物学靶点和分子机制，为结直肠癌的诊断和治疗提供了新的建议和理论支持。

Abstract: Objective: To explore the potential therapeutic target and molecular mechanism of berberine in the treatment of colorectal cancer. Methods: Swiss Target Prediction, Comparative Toxicogenomics, TargetNet and Binding databases were used to screen the active components and gene targets of berberine. In parallel, WGCNA analysis was performed based on the GSE90627 dataset from the GEO database to identify relevant targets for colorectal cancer. The disease-drug gene set is ob-tained by drawing a Venn diagram. A protein interaction network (PPI) was constructed using STRING database. Topology analysis was performed using Cytoscape software to obtain key cluster. The key cluster genes were analyzed for functional enrichment. Finally, virtual validation was per-formed through an online molecular docking website. Results: After a series of bioinformatics anal-yses, key clusters containing 35 genes were obtained. Enrichment analysis showed that key cluster genes were mainly enriched in p53 and FOXO signal pathways. Four of the most important biologi-cal targets were identified: WD repeat protein 74 (WDR74), cyclin dependent kinase 4 (CDK4), au-rora kinase A (AURKA), and cyclin E1 (CCNE 1). Molecular docking results showed that all four tar-gets could bind to berberine stably. Conclusion: In this study, biological targets and molecular mechanisms of berberine in the treatment of colorectal cancer were screened and identified by bi-oinformatics and network pharmacology methods, which provided new suggestions and theoretical support for the diagnosis and treatment of colorectal cancer.

1. 引言

结直肠癌(colorectal cancer, CRC)是一种常见的消化系统恶性肿瘤，在最近的全球癌症发病率和死亡率排名中排名居高不下 [1] 。它的高发病率和致死率对人类生命和福祉构成严重威胁 [2] 。结直肠癌复杂的发病机制和发生发展涉及大量与肿瘤相关的基因和信号通路的改变。此外，缺乏有效的临床筛查生物标志物使得早期诊断变得困难。据研究，近20%的结直肠癌患者在确诊时已发生不同程度的转移，5年总生存率不足15% [3] [4] 。放疗、化疗、手术和免疫治疗是目前结直肠癌的主要治疗手段 [5] [6] 。然而，由于它们存在一些局限性，如全身毒性、药物抗性、低选择性、低效性和额外的副作用等等，对人类生命、健康和经济财产构成了重大威胁 [7] 。因此，寻找更安全、成本更低和更有效的结直肠癌治疗方法成为了普遍趋势。

近年来，作为一种古老的东方医学，中医药在全球范围内作为一种辅助医学得到了广泛认可。目前，使用具有抗癌作用的中草药提取物进行化疗，是除手术和放疗之后的第三大癌症常用治疗方法 [8] [9] [10] 。中医对结直肠癌的认识可以追溯到几千年前。目前已研制出许多治疗结直肠癌疗效确切的中药汤剂，包括先练解毒汤 [11] 、左金丸 [12] 、扶正解毒汤 [13] 等，它们均含有相同的中药黄连，其主要药理成分为小檗碱。先前的研究表明小檗碱对结直肠癌疗效确切 [14] [15] [16] ，然而，这些效果的具体机制仍然不明确。

基于生物信息学的网络药理学是一种在生物分子网络中绘制药物靶点和疾病相关分子的新策略 [17] 。随着生物信息学和药理学的快速发展，网络药理学研究可视化了药物、靶点和疾病之间的相互关系 [18] 。本研究基于加权基因共表达网络分析(WGCNA)，探讨小檗碱治疗结直肠癌的治疗靶点、分子机制和所涉及的相关通路。

2. 材料与方法

2.1. 数据收集和预处理

在本研究中，我们从PubChem数据库(https://pubchem.ncbi.nlm.nih.gov)中获得了小檗碱的化学结构式和标准SMILES。然后，我们使用获得的小檗碱的化学结构式在TargetNet数据库(http://www.example.com)、Binding数据库(https://www.bindingdb.org/)、Comparative Toxicogenomics数据库(http://ctdbase.org/)和SwissTargetPrediction数据库(http://www.swisstargetprediction.ch/)中搜索和整理小檗碱的潜在药理靶点。最后，使用R软件(4.2.1版本)的“venn”程序包对其进行可视化。结直肠癌的表达谱数据集GSE 90627下载自基因表达综合(GEO)数据库(https://www.ncbi.nlm.nih.gov/geo/)。该数据集包含32个结直肠癌活检样本和96个正常组织样本，每个样本在液氮中速冻并石蜡包埋以便于后续分析。所有组织样本均经组织学证实，并获得临床病理学信息。

2.2. 数据标准化和差异表达基因的识别

使用“limma”程序包的“normalizeBetweenArrays”函数进行标准化 [19] 。使用“limma”软件包对标准化数据集进行样品之间的差异基因分析，以调整后的p < 0.05和|log2FC| > 1作为显著性筛选标准。然后使用“pheatmap”和“gglot2”包绘制热图和火山图，以可视化结果 [20] 。

2.3. 加权基因共表达网络分析(Weighted Gene Co-Expression Network Analysis, WGCNA)的构建及模块鉴定

为了探索基因之间潜在的调控关联，我们使用R软件的“WGCNA”包构建基因共表达网络。首先，为保证网络结果的可靠性，以平均片段数(FPKM) > 0.5作为筛选标准，对正常值样本进行筛选。随后，使用“flash Clust”工具包进行聚类分析，保留在指定阈值下最适合聚类的样本。然后，利用“pick Soft Threshold”函数选择最优加权系数β值的软阈值，建立无标度网络，从而将相似度矩阵转化为邻接矩阵。将邻接矩阵转换为拓扑重叠矩阵，并计算相应的不相似度。此外，将模块裁剪高度设为0.3，最小模块数设为100，然后采用动态树裁剪法从层次聚类树中识别模块 [21] 。最后筛选出与结直肠癌关系最密切的模块，并将该模块中的基因与差异分析筛选出的基因进行交集，作为疾病相关基因。

2.4. 蛋白互作网络的建立和关键簇的识别

我们对疾病和药物靶点取交集，然后将其上传到STRING数据库(https://www.string-db.org/)以构建蛋白互作网络(PPI)。通过Cytoscape软件(3.9.1版)的分子复合物检测(MCODE)鉴定PPI网络的关键簇 [22] 。对正常组和肿瘤组之间筛选的关键簇中的基因表达进行Wilcoxon秩和检验。所有检验均为双侧检验，认为p值小于0.05具有统计学显著性。

2.5. 关键簇的功能丰富分析

为探讨小檗碱治疗结直肠癌的生物学机制，我们对关键基因簇进行了富集分析。我们分别进行了基因本体论(GO)和京都基因与基因组百科全书(KEGG)富集分析。探讨其潜在的生物学过程、功能靶点、信号通路，以及两者之间的潜在关系。所有分析均通过R软件的“clusterProfiler”软件包进行，调整后的q值 < 0.05被设定为显著富集阈值。

2.6. 分子对接

为了验证小檗碱与预测的核心靶点的结合能力，从PubChem数据库检索小檗碱的MOL2分子结构，从RCSB数据库(https://www.rcsb.org/)检索核心靶点的三维晶体结构 [23] ，然后使用分子一键对接网站(https://mcule.com/apps/1-click-docking/)进行受体和配体的在线分子对接。

3. 结果

3.1. 数据收集和处理

我们从4个数据库中共收集了总共190个小檗碱的潜在药理靶点，其结果在维恩图中呈现(图1A)。此外，完成归一化后的GSE90627数据集所有样本居中，数值分布合理，表明该微阵列数据集质量较高，交叉可比性良好。

3.2. 差异基因的识别

使用调整后的p < 0.05和|log2 FC| > 1作为筛选条件，在GSE90627数据集中共识别出1976个差异基因，其中包含852个上调基因和1124个下调基因(图1B)。此外，分别选取排名前30的上、下调基因绘制热图(图1C)。

Figure 1. Potential pharmacological targets of berberine and differential genes in colorectal cancer

图1. 小檗碱的潜在药理靶点与结直肠癌差异基因

3.3. 加权共表达网络分析

首先，我们通过Pearson相关分析筛选了总共16,288个用于WGCNA分析的基因。然后，根据“WGCNA”包的“pick Soft Threshold”函数的计算结果，我们选择5作为最佳软阈值。此外，为了使共表达网络符合无标度原理，本研究选取β = 7 (R² = 0.9)构建基因共表达网络(图2A、B)。最后，我们将相似度矩阵转换为拓扑重叠矩阵。接下来，在拓扑重叠矩阵的基础上，我们使用前面提到的标准来合并相似度较高的模块，最终得到总共17个模块(图2C)。最后，根据模块–性状关联分析结果，我们选择与结直肠癌关联最强的黑色模块(总共885个基因)，并与先前获得的差异基因(总共1976个基因)取交集，最终获得169个疾病基因(图2D)。

Figure 2. Weighted coexpression network analysis

图2. 加权共表达网络分析

3.4. 关键簇的鉴定和功能富集分析

使用STRING在线数据库，我们使用疾病–药物基因集构建了一个由336个节点和2020条边组成的蛋白质–蛋白质相互作用网络，并使用Cytoscape对其进行可视化。然后，使用MCODE插件识别基因簇模块，并选择得分最高(得分：15.167，37个节点和546条边)的模块作为关键模块(图3A)。对疾病–药物基因集的关键簇进行进一步的功能富集分析。GO富集分析结果显示，本体上共有424个基因被显著富集(q < 0.05)，其中生物过程(BP)基因353个，细胞组分(CC)基因49个，分子功能(MF)基因22个。如图3B所示，它们主要富集在ncRNA加工、转移酶复合物和蛋白丝氨酸激酶活性等过程中。KEGG富集分析结果显示，共有41条途径被显著富集(q值 < 0.05)，包括细胞衰老、p53信号通路和FOXO信号通路等等(图3C)。

Figure 3. Disease-pharmacogene key clusters and GO analysis results

图3. 疾病–药物基因关键簇和GO分析结果

3.5. 分子对接验证

我们利用筛选的最后4个核心靶点与小檗碱(MOL 001454)进行分子对接验证。结果显示，4个核心靶点与小檗碱的结合能得分均大于7.0 kcal·mol⁻¹。具体而言，WDR74 (PDB ID: 5SUI)与小檗碱的结合能评分为−8.8 kcal·mol⁻¹ (图4A)，CDK 4 (PDB ID: 7 OXW)与小檗碱的结合能评分为−7.3 kcal·mol⁻¹ (图4B)，AURKA (PDB ID: 6VPM)与小檗碱的结合能评分为−8.4 kcal·mol⁻¹ (图4C)，CCNE 1 (PDB ID: 6 GUK)与小檗碱的结合能评分为−9.5 kcal·mol⁻¹ (图4D)。这些发现证明小檗碱实际上能够稳定地结合所评估的四个主要靶点，以及小檗碱能够有效地靶向结直肠癌治疗中的WDR74、CDK4、AURKA和CCNE1。

Figure 4. Verification of molecular docking between the core target and berberine (MOL 001454)

图4. 核心靶点与小檗碱(MOL 001454)分子对接验证

4. 讨论

近年来，发展中国家的结直肠癌发病率和死亡率正逐年上升。小檗碱是毛茛科植物黄连的干燥根茎，是许多常用的治疗肠道疾病的中药处方中的主要成分。本研究利用网络药理学和生物信息学方法，对黄连中小檗碱与结直肠癌的关系进行了研究，以期为进一步了解结直肠癌的发病机制和寻找新的治疗方法提供理论依据。在本研究中，我们首先对GEO数据库的GSE90627数据集进行了差异分析和WGCNA分析。在获得与结直肠癌关联性最强的差异基因和黑色模块后，将差异基因与黑色模块中的基因取交集，获得169个疾病基因。随后，我们将疾病基因与小檗碱的所有潜在靶点进行了比较，这些靶点先前已在四个不同的数据库中进行了搜索，以进行联合分析。我们为获得的疾病–药物基因组构建了PPI网络，并确定了最关键的基因簇(得分：15.167)。结果显示，小檗碱对大肠癌的治疗作用可能通过多种途径实现，其中p53和FOXO信号通路的表达量最高。先前的研究表明，p53信号通路主要负责调节TP 53的表达，TP 53通过上调DNA修复蛋白、保持基因组稳定性和诱导细胞凋亡来控制癌症的进展 [24] [25] 。此外，TP 53还可以激活miR-22的表达，miR-22已被证明显著抑制AP-1结合DNA的能力，并且是一种关键的肿瘤抑制miRNA [25] [26] 。以前的研究已经证明，功能失调的FOXO蛋白是肿瘤发生和发展的重要因素 [27] [28] 。值得注意的是，FOXO途径经常与其他代谢途径结合，以抑制或阻断癌症的进展，作为作用于肿瘤的手段。通过调节FOXO及其相关通路，许多药物，如紫杉醇、拉帕替尼、表柔比星、顺铂和多柔比星，在临床上用于治疗各种癌症 [26] [29] 。提示小檗碱可能通过上述两种机制作用于直肠癌组织。此外，根据对小檗碱治疗结直肠癌的靶点的GO富集分析结果，小檗碱的抗CRC影响可能与rRNA加工和代谢、ncRNA加工和细胞有丝分裂进程有关。

随后，我们筛选出与小檗碱抗结直肠癌最相关的4个标志物：WDR 74、CDK 4、AURKA和CCNE 1。通过对临床样本表达的差异分析和分子对接验证，我们认为上述4个标志物可能是小檗碱治疗结直肠癌的关键靶点。WD重复蛋白74 (WDR 74)是一种核糖体装配因子，主要参与rRNA合成和核糖体生物发生 [30] 。Li等发现，WDR 74通过调控RPL 5-MDM 2-p53通路刺激p53网络下游分子，在功能上促进黑色素瘤的发生和转移 [31] ，这与我们之前KEGG富集分析的结论一致。此外，许多先前的研究已经得出结论，WDR 74是通过调节连环蛋白信号通路而导致结直肠癌发生和进展的直接原因，这与我们的发现一致 [32] 。细胞周期蛋白依赖性激酶(CDK) 4有助于肿瘤的发生。目前，临床上使用的许多抗肿瘤药物通过促进CDK 4降解来靶向恶性肿瘤 [33] 。Assi教授等人认为CDK 4与脂肪瘤的病因学密切相关，并指出CDK 4抑制剂是脂肪肉瘤最合理的治疗靶点 [34] 。Aurora激酶A是Aurora激酶家族的三个成员之一，是一种有丝分裂丝氨酸/苏氨酸激酶，是一种重要的有丝分裂调节因子 [35] 。多项研究报告AURKA的异常激活或过表达是结直肠癌发生或进展的主要原因 [36] [37] 。此外，有研究报告了某些草药化合物成分通过p53信号通路治疗CRC的机制，这与我们的发现一致 [38] 。目前尚未有报道CCNE 1与结直肠癌的相关性，但已在其他癌症中对其进行了研究 [39] 。我们未来的研究还将探讨CCNE 1与结直肠癌的药理学关系。此外，这项研究同样存在一些局限性。首先，在数据收集方面，为了补充疾病靶点并结合患者结肠粘膜样本验证差异基因，我们可以进一步搜索其他疾病数据库。其次，我们还需要进行相关细胞和动物实验以验证相关的靶点和信号通路。