青岛大学附属医院老年医学科，山东 青岛

收稿日期：2023年11月11日；录用日期：2023年12月4日；发布日期：2023年12月13日

摘要

目的：通过生信分析筛选原发性双侧大结节肾上腺增生的核心基因，探索疾病治疗的新靶点。方法：本研究从基因表达数据库(Gene Expression Omnibus, GEO) (http://www.ncbi.nlm.nih.gov/geo)检索并下载了与原发性双侧大结节肾上腺增生相关的转录组测序数据和表达谱数据集GSE171558。我们筛选了差异表达基因(Differentially Expressed Genes, DEGs)，并对该数据集进行了加权基因共表达网络分析(Weighted Gene Co-Expression Network Analysis, WGCNA)。对蓝色模块基因进行了基因本体论(Gene Ontology, GO)、京都基因和基因组百科全书(Kyoto Encyclopedia of Genes and Genomes, KEGG)通路富集分析和基因集富集分析(Gene Set Enrichment Analysis, GSEA)。根据差异表达基因进行了蛋白质–蛋白质相互作用网络(Protein-Protein Interaction Network, PPI)。我们选择了蓝色模块中的中枢基因和PPI中的枢纽基因的重叠基因作为PBMAH的最终中枢基因。并且我们在GSE25031数据集中验证了最终枢纽基因。结果：蓝色基因模块(677个基因)主要富集在脂质代谢领域，与PBMAH具有最高的相关系数。通过差异分析，我们筛选出487个DEGs，其中包括231个上调基因和256个下调基因。蛋白质–蛋白质相互作用网络分析鉴定出30个中枢基因。GPC4和VCAN被确定为PBMAH的最终中枢基因。与正常对照组相比，GPC4在PBMAH组中的表达显著下调，而VCAN与正常组相比显著上调。GSEA数据分析显示，VCAN与PI3K-Akt信号通路、磷脂酶D信号通路、Rap1信号通路、Ras信号通路、MAPK信号通路等相关联。GPC4与癌症相关通路、Rap1信号通路、PI3K-Akt信号通路、Wnt信号通路等相关。GSE25031的原始数据验证了分析结果。结论：基于生物信息学分析初步筛选出GPC4和VCAN可能参与PBMAH的致病和发展。为进一步研究原发性双侧大结节肾上腺增生的诊断及治疗提供了新的方向。

关键词

原发性双侧大结节肾上腺增生，库欣综合征，WGCNA，枢纽基因，GEO，GSEA，PPI，KEGG，GO， VCAN

Identifying and Validating GPC4 and VCAN as Hub Genes in Primary Bilateral Macronodular Adrenal Hyperplasia by WGCNA and DEG

Yinfei Xu, Yuhe Jiang, Hui Yan, Wenjing Feng, Xiaotong Pan, Caixia Cao^*

Department of Geriatrics Medicine, Affiliated Hospital of Qingdao University, Qingdao Shandong

Received: Nov. 11^th, 2023; accepted: Dec. 4^th, 2023; published: Dec. 13^th, 2023

ABSTRACT

Objective: To identify hub genes in primary bilateral macronodular adrenal hyperplasia (PBMAH) through bioinformatics analysis and explore novel targets for disease treatment. Methods: This study searched and downloaded the transcriptome sequencing data and expression profile dataset GSE171558 related to primary bilateral adrenal macronodular hyperplasia from the gene expression omnibus, GEO, http://www.ncbi.nlm.nih.gov/geo). We filtered the differentially expressed genes (DEGs) and performed weighted gene coexpression network analysis (WGCNA) on this dataset. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment Analysis and Gene Set Enrichment Analysis (GSEA) were performed for the blue module genes. Protein-protein interaction network (PPI) analysis was performed based on the differentially expressed gene. We selected the overlapping genes of the hub gene in the blue gene module and the hub gene in PPI as the final hub gene of PBMAH. And we verified the final hub gene in the GSE25031 dataset. Results: The blue gene model (677 genes) is mainly enriched in lipid metabolism, with the highest correlation coefficient with PBMAH. Through differential analysis, we screened out 487 DEGs, including 231 up-regulated genes and 256 down-regulated genes. PPI analysis identified 30 hub genes. GPC4 and VCAN were identified as the final hub genes of PBMAH. The raw data of GSE25031 verified the analysis results. The expression of GPC4 was significantly down-regulated in the PBMAH group compared with the normal control group, and VCAN was up-regulated considerably compared with the normal group. Analysis of GSEA data showed that VCAN was connected to PI3K-Akt signalling pathway, Phospholipase D signalling pathway, Rap1 signalling route, Ras signalling pathway, MAPK signalling pathway, etc. GPC4 was associated to cancer-related Pathways, Rap1 signalling pathway, PI3K-Akt signalling pathway, Wnt signalling pathway, etc. Conclusions: Based on bioinformatics analysis, GPC4 and VCAN have been initially identified as potential players in the pathogenesis and development of PBMAH. This provides a new direction for further research on the diagnosis and treatment of primary bilateral macronodular adrenal hyperplasia.

Keywords:Primary Bilateral Macronodular Adrenal Hyperplasia, Cushing Syndrome, WGCNA, Hube Genes, GEO, GSEA, PPI, KEGG, GO, VCAN

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

库欣综合症(Cushing syndrome, CS)通常是由于肾上腺过度分泌皮质醇，从而抑制来自垂体前叶的促肾上腺皮质激素(adrenocorticotropic hormone, ACTH)的释放，导致血浆ACTH水平降低 [1] 。可根据结节的大小分为原发性双侧大结节肾上腺增生(primary bilateral macronodular adrenal hyperplasia, PBMAH)和微小双侧肾上腺皮质增生(micronodular bilateral adrenocortical hyperplasia, MiBAH) [2] 。PBMAH患者体内过量的皮质醇通常是由于激活肾上腺皮质异常G蛋白偶联受体(G protein-coupled receptors, GPCR)而导致肾上腺细胞增生 [3] 。当异位受体如胃抑制性肽受体(Gastric Inhibitory Peptide Receptor, GIPR)、β-肾上腺素受体(Beta-adrenergic Receptor, β-AR)与配体结合时，将激活细胞内信号转导途径，促进皮质醇的分泌 [4] 。近年来，ARMC5基因也通过综合全基因组方法被确认为表观上散发或家族性PBMAH的常见致病基因之一 [5] ，胚系ARMC5突变是常见的遗传缺陷 [6] 。ARMC5的致病性变体导致了20%~25%的PBMAH,通常会导致更严重的表型 [3] [5] [7] 。Linjiang Lao [8] 等人认为ARMC5突变会导致RNA聚合酶II (RNA Polymerase II, Pol II)池的扩大，而导致效应基因子集失调，这种扩大的Pol II池和基因失调与肾上腺增生相关。ARMC5被发现是PBMAH中频繁发生突变的抑癌基因，它的失活会通过增加转录调控因子NRF1的表达，能提高清除活性氧自由基的因子表达，进而调节肾上平衡腺皮质细胞中氧化还原 [9] 。另有研究揭示了ARMC5还控制类固醇激素的合成和细胞存活 [9] [10] 。WGCNA被认为是一种处理基因表达数据的高效方法，有助于发现与临床病例特征相关的基因，现广泛应用于各种癌症的遗传分析中 [11] [12] [13] 。我们将在本研究中使用WGCNA来识别潜在的治疗和干预PBMAH的靶点基因。

2. 材料及方法

2.1. 数据析取

来自GEO公共数据库的GSE171558数据集，包括三名受影响的兄弟的两个正常肾上腺和两个结节，使用WGCNA和DEG筛选分析该数据集，以确定与PBMAH相关的关键模块和枢纽基因。该数据集基于芯片平台GPL6244 Affymetrix 1.0 ST Array。

2.2. 患者基本资料

III-1、III-2和III-3是一个家族的三个男性兄弟。III-1是由于PBMAH导致晚期库欣综合症，已行双侧肾上腺切除手术。他的两个兄弟，III-2和III-3，患有轻度高皮质醇症，已行单侧肾上腺切除手术。他们的年龄、性别和诊断情况如下所示(表1)。本研究流程图见图1。

图1. 本研究流程图

表1. 患者基本资料

注：PBMAH：原发性双侧大结节性肾上腺增生。

2.3. 加权基因共表达分析

使用R软件“WGCNA”包执行了加权基因共表达分析 [14] 。基因的模块关系值(Module membership, MM)体现了基因与模块的相关程度。利用 WGCNA包的signed KME计算模块关系值，选取目标模块中 MM的绝对值排名前5%的转录因子作为该模块的核心基因，并筛选与关键转录因子连通性(weight值)前50的关联节点进行基因调控网络构建，将数据导入Cytoscape软件 [15] ，构建基因共表达网络图。我们计算了各组之间每个基因差异表达的显著性检验P值，并以10为底取对数作为基因显著性指数(gene significance, GS)。根据GS和MM值进行筛选。MM通过R软件的WGCNA包中的biserial.cor函数计算，代表基因与特征之间的关联程度。我们筛选了GS和MM值都大于0.9的模块基因作为中枢基因。差异表达基因根据适当的β值构建成邻接矩阵。随后，我们将其转化为拓扑转移矩阵(topological transition matrix, TOM)以识别基因模块。最后，我们选择与PBMAH产生高度关联的基因进行后续分析。为避免性别和年龄对基因表达的影响，我们分析了样本的基因与性别及年龄之间的相关性。

2.4. 富集分析

我们使用基因本体论富集分析(Gene Ontology, GO)和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes, KEGG)通路富集分析。由R语言中“clusterProfiler”包 [16] 完成，筛选条件是P < 0.05。基因集富集分析(Gene Set Enrichment Analysis, GSEA)用于关键模块基因的富集分析。根据中枢基因表达的中位数，将样本分为高表达组和低表达组，并按差异倍数值从大到小排序。使用R包“clusterProfiler”对基因进行功能注释，以全面探讨这些差异基因的功能相关性。

2.5. 差异表达基因的筛选

使用R软件包 [17] 中的“DESeq2”包对数据集进行差异分析，差异基因筛选条件为|logFC| > 1 &adj.P.Val < 0.05。使用“Ggplot2”可视化差异基因。筛选后，使用STRING软件在StringData 11.5 (https://string-db.org/cg/input.pl)中构建差异表达基因的蛋白质-蛋白质相互作用网络(Protein-protein interaction network, PPI)。设置最小的互动得分为0.9。使用Cytoscape软件来可视化PPI网络 [18] 。cytoHubba软件中的分子复杂检测(Molecular Complex Detection, MCODE)被用来筛选出前30个基因作为关键基因 [19] 。

2.6. 枢纽基因的验证

在PPI网络和WGCNA中均存在的基因被选择为最终中枢基因。我们设定以下条件来找到GEO网站上合适的数据集来验证枢纽基因：(1) 具有PBMAH患者肾上腺腺体的测序数据；(2) 具有完整的基因矩阵信息；(3) 至少三个样本。GEO公共数据库下载的GSE25031数据集符合这些标准，包括正常对照组和七个PBMAH患者的肾上腺结节数据。

3. 结果

3.1. 加权基因共表达分析

当相关系数为0.8时，选择β = 12作为软阈值(R2为0.99) (图2(B)和图2(C))，通过样本聚类来查看异质性(图2(A))。WGCNA分析生成了35个不同灰度颜色显示的基因模块(图3(A))。基于模块基因与特征之间的相关性分析结果制作相关性热图(图3(B)和图3(D))。结果提示这些基因模块均与PBMAH的发生和发展有关。蓝色基因模块(677个基因)与PBMAH具有最高的相关系数(CC为0.88，P = 4.70E−05)，呈正相关。随机选择了400个基因进行热图可视化(图3(C))。相关性分析结果显示，蓝色模块基因与年龄无关(P > 0.05) (图3(B))。蓝色基因模块被用作后续分析的关键基因模块。

图2. 数据验证和过滤器的最佳软阈值。(A) 通过样本聚类来查看异质性。(B) 软阈值非尺度拟合指标分析。(C) 软阈值平均连通性分析

图3. 基于GSE171558数据集构建的共表达网络。(A) 基因树状图。(B) 模块–特征关系的热图。(C) 热图可视化了在网络中随机选择的400个基因来描述TOM。(D) 特征基因树状图和热图的结合

3.2. 蓝色模块基因的功能分析

我们对该模块中的基因进行了GO富集分析和KEGG富集分析(图4(A)，图4(B))。在生物过程方面，GO分析显示这些基因主要富集在这些细胞生物过程(biological processes, BP)，包括丙酮酸合成乙酰辅酶A，乙酰辅酶A生物合成与加工和三羧酸循环等方面。细胞组分(Cellular components, CC)包括线粒体基质、氧化还原酶复合物和二氢脂酰亚硫酸脱氢酶复合物。分子功能(Molecular functions, MF)包括作用于供体的醛或氧羰基，以NAD或NADP为受体的氧化还原酶活性和作用于供体的醛或氧羰基和铁硫簇结合的氧化还原酶活性。KEGG分析确定柠檬酸循环、丙酮酸代谢和糖酵解/糖异生途径明显富集。这证实了这些基因在散发性PBMAH的发展中起着关键作用。GSEA数据分析显示VCAN与PI3K-Akt信号通路、黏附斑、黑色素瘤、缝隙连接、前列腺癌、磷脂酶D信号途径、Rap1信号途径、Ras信号途径、MAPK信号途径和癌症中胆碱代谢等有关(图4(C))。GSEA数据分析显示GPC4与轴突引导、ECM受体相互作用、癌症途径、Hedgehog信号通路、黏附斑、基底细胞癌、Rap1信号通路、调控肌动蛋白细胞骨架、PI3K-Akt信号通路和Wnt信号通路等有关(图4(C))。

图4. 枢纽基因的富集分析。(A) 功能聚类-GOKEGG (联合logFC)-圈图。内圈的每个柱子对应一个条目，高度为p.adj的相对大小。柱子对应填充的颜色代表条目对应的z-score值。z-score = (Up − Down)/ $\sqrt{\text{Counts}}$ ，Counts代表条目对应的分子总数。z-score为正，说明对应的条目可能是正调节，如果为负，对应条目可能是负调节。(B) 功能聚类-GOKEGG (联合logFC)-气泡图，圆点直径越大表示富集的基因数越多。(C) VCAN和GPC4的基因集富集分析

3.3. 差异表达基因的筛选

在使用GSE171558数据集进行筛选后，共发现了487个差异表达基因(DEGs)，其中231个基因表达上调，256个基因表达下调(图5(A)和图5(B))。当k分数大于2时，使用MCODE软件包筛选出了共30个核心基因，分布在三个簇内。这三个簇分别包含45个节点(图5(C))、44个节点(图5(D))和19个节点(图5(E))。

图5. 差异基因的筛选(A, B)通过筛选GSE171558数据集，得到的火山和热图。(C, D, E)使用MCODE插件进行基因簇过滤

3.4. 枢纽基因的验证

MYH11、CNN1、MYOCD、SMTN、NOTCH1、VCAN和GPC4，同时存在于WGCNA分析和PPI分析结果中，因此被视为最终的核心基因。为了验证结果，使用GSE25031数据集 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse25031)对这七个基因进行了验证，结果显示GPC4显著下调，VCAN显著上调(图6)。

图6. 外部数据集的验证(A, B) GPC4、VCAN在GSE25031数据集中表达量箱线图

4. 讨论

以往的研究表明体细胞基因突变参与了PBMAH的发病机制，如腺瘤性息肉病基因(APC) [20] 、多发性内分泌肿瘤类型1 (MEN1) [21] 和富马酸酶(FH) [22] 。最近，通过全基因组方法识别的包含ARMC5的突变也被认为是散发或家族性PBMAH的常见可能的致病基因 [23] [24] [25] 。PBMAH目前病因尚不明确，基因测序技术较常被运用于分析PBMAH的致病基因，而WGCNA和其他生物信息学分析方法的联合运用尚未被用于筛选PBMAH的枢纽基因 [26] 。我们使用WGCNA方法分析PBMAH患者肾上腺的生物信息，以确定PBMAH的枢纽基因。为了消除混杂因素，使用相关性分析来避免性别和年龄对基因表达的影响。结果显示蓝色模块基因与患者是否患PBMAH有关，但与年龄无关(P > 0.05)。因此，本研究可以避免年龄对基因表达量的影响。本研究根据WGCNA和MCODE分析确定的最终核心基因GPC4和VCAN被定义为关键的核心基因。其中，富集分析结果主要聚集在脂质代谢，表明异常的脂质代谢可能是导致PBMAH患者皮质醇异常产生的病因。GPC-4与细胞膜相关联的肝素硫酸糖蛋白 [27] ，与胰腺癌细胞的迁移、增殖、分化和形态生成有关 [28] 。其表达下调可能导致乳腺癌细胞的增加迁移和增殖 [29] 。GSEA分析结果显示，GPC4与癌症相关的途径、Rap1信号通路、PI3K-Akt信号通路和Wnt信号通路等有关。Rap1是一种调控代谢的蛋白质，包括细胞膜受体的信号传导、细胞分裂所需的细胞骨架重排等 [30] ，已被证明与肿瘤形成有关 [31] 。Rap1信号通路的激活可以同时导致炎症介质的生成，这对于调节细胞外基质和影响纤维化的过程至关重要 [32] 。我们猜测Rap1表达的异常可能影响肾上腺周围组织的炎症反应以及肾上腺细胞的增殖和纤维化 [33] 。根据先前的研究，LPS可以激活p38 MAPK和PI3K/Akt信号通路，刺激小鼠肾上腺皮质中的Y1细胞合成类固醇 [34] 。本研究GSEA分析显示，GPC4与PI3K/Akt信号通路相关，由此我们猜测人体中可能也存在此途径介导的类固醇异常合成。本研究显示，VCAN在PBMAH患者中表达上调。既往已经证实VCAN-AS1的过表达可以加速细胞的增殖、迁移、侵袭和肿瘤生长，同时抑制细胞凋亡 [35] 。而本研究GSEA数据分析显示，VCAN与PI3K-Akt、磷脂酶D、Rap1、Ras、MAPK信号通路相关联。研究发现，磷脂酶D对AngII引起的肾上腺类固醇反应和皮质醇释放至关重要 [36] 。此外，皮质醇还可能通过激活PI3K-Akt信号通路来增加牛子宫内上皮细胞(Bovine Endometrial Epithelial Cells, BEECs)的增殖 [37] 。PI3K/Akt被证实能增加长期缺氧绵羊胎儿肾上腺皮质细胞中皮质醇生物合成 [38] 。因此，我们推测VCAN可能通过上述途径刺激肾上腺增殖及调节皮质醇生成 [39] 。Rap1已被证明可通过刺激PI3K/AKT信号，增强肝细胞癌(Hepatocellular Carcinoma, HCC)的生长和转移 [40] ，并且可能影响肾上腺细胞增殖 [41] 。以往研究已证明类固醇急性调节蛋白可通过MAPK途径参与类固醇合成 [42] 。提示我们PBMAH患者MAPK途径的异常可能参与了类固醇异常生成 [43] 。综上所述，我们的研究通过WGCNA和DEG分析确定了GPC4、VCAN是PBMAH的最终枢纽基因。但是，这些枢纽基因关于PBMAH发病的具体机制仍需要进一步研究来证明。本研究为进一步研究PBMAH的发病机制及治疗靶点提供了新的方向。

文章引用

徐音飞,江予赫,闫 慧,冯文静,潘晓彤,曹彩霞. 通过WGCNA和DEG确定并验证GPC4和VCAN作为原发性双侧大结节肾上腺增生的枢纽基因
Identifying and Validating GPC4 and VCAN as Hub Genes in Primary Bilateral Macronodular Adrenal Hyperplasia by WGCNA and DEG[J]. 临床医学进展, 2023, 13(12): 19010-19022. https://doi.org/10.12677/ACM.2023.13122674

参考文献