基于WGCNA的尿道下裂相关基因鉴定分析

刘晓雅 常孟孟 马文岳 高洪杰 孙丰银

[摘要] 目的 利用基因共表達权重网络（weighted gene co-expression network analysis，WGCNA）探索与尿道下裂发病相关的潜在基因。方法 利用WGCNA算法将数据集GSE35034处理并构建基因共表达权重网络，筛选出与尿道下裂相关性最高的模块，通过基因本体论（gene ontology，GO）、京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）对模块中的基因进行富集分析。同时进行差异分析，筛选差异基因。将差异基因导入String数据库，利用Cytoscape软件，找出网络中枢纽基因。将以上3 种方法的结果系统分析，筛选出交集中的核心基因。利用外部数据集对核心基因进行mRNA表达变化及受试者操作特征（receiver operating characteristic，ROC）曲线诊断验证。结果 基于WGCNA方法得到15个共表达模块。通过差异分析获得了93个满足条件的共同差异基因。通过Cytoscape软件分析最终得到10个核心基因。最后3种方法得到2个交集基因：FBXL16、SYNDIG1。ROC曲线验证了交集基因在尿道下裂患者与正常人中表达存在差异。结论 本研究通过WGCNA获得了2个与尿道下裂具有显著相关性的关键基因，可用于尿道下裂发病及治疗的探索。

[关键词] 尿道下裂；基因共表达权重网络；差异分析；富集分析

[中图分类号] R726.9      [文献标识码] A      [DOI] 10.3969/j.issn.1673-9701.2024.14.002

Analysis of key gene related to hypospadias based on gene co-expression weight network analysis

LIU Xiaoya1， CHANG Mengmeng1， MA Wenyue1， GAO Hongjie2， SUN Fengyin1，3

1.Department of Pediatric Surgery， Qilu Hospital of Shandong University， Jinan 250012， Shandong， China; 2.Department of Pediatrics， Qilu Hospital of Shandong University， Jinan 250012， Shandong， China; 3.Center for Post-Doctoral Studies of Shandong Qidu Pharmaceutical Co.，LTD， Jinan 250012， Shandong， China

[Abstract] Objective To explore potential genes associated with the pathogenesis of hypospadias using weighted Gene co-expression network analysis （WGCNA）. Methods The WGCNA algorithm was used to process the hypospadias-related dataset GSE35034， and then a gene co-expression weight network was constructed to screen the modules with the highest correlation with hypospadias， and the genes in the modules were enriched and detected by gene ontology （GO）， Kyoto Encyclopedia of Genes and Genomes （KEGG）. Differential analysis was also performed to screen out differential genes. The differential genes were imported into the String database. Using Cytoscape software， the hub genes in the network were identified. The results screened by the above three methods were combined and analyzed， and the core genes in the intersection set were screened. Using the external dataset GSE121712， the core genes were verified by mRNA expression changes and subject work characterization receiver operating characteristic （ROC） curve diagnosis. Results Fifteen co-expression modules were obtained based on the WGCNA method. Ninety-three common differential genes  meeting the conditions were obtained by differential analysis. Ten core genes were finally obtained by Cytoscape software analysis. Finally MEbrown module， differential genes and the 10 core genes yielded a total of 2 intersecting genes： FBXL16 and SYNDIG1. ROC curves verified that the intersecting genes were differentially expressed in patients with hypospadias versus normal subjects . Conclusion In this study， two key genes with significant correlation with hypospadias were obtained by WGCNA， which may be used for the early diagnosis and treatment of hypospadias after further study.

[Key words] Hypospadias; Weighted gene co-expression network analysis; Differential analysis; Enrichment analysis

尿道下裂是最常见的男性泌尿生殖系统畸形之一，但其发病机制仍不明确。近年来，随着外科手术技术不断进步，已出现多种手术方式，但治疗效果仍不满意；近几十年来发病率有上升趋势，因此人们对尿道下裂的发病原因也更加关注[1]。尿道下裂的遗传易感性具有公认的家族性聚集性，患有尿道下裂的男孩的男性亲属更容易患这种疾病。当遗传易感性与抗雄激素药物暴露相结合时，基因和环境风险因素可能相互作用超过阈值，导致这种出生缺陷的发生[2]。因此，尿道下裂发育畸形和生殖疾病的发生，存在基因和环境的共同作用，并且携带遗传易感基因合并环境内分泌干扰物的暴露可能是发病的主要原因。

在基因组时代，随着高通量测序技术发展和基因表达数据库的完善，为尿道下裂的病因学研究提供了新的途径。基因表达数据库（gene expression omnibus，GEO）是当今最全面的基因表达数据库，目前储存了高达208万样本的表达谱数据。本研究从GEO上选择了包含3种人类包皮成纤维细胞基因表达谱的芯片，进行加权基因共表达网络分析。

WGCNA是一种基于基因表达数据集检测，研究高度相关基因簇的常用方法，通过构建基因表达网络，筛选生物标志物或治疗靶点[3-4]。WGCNA可根据基因表达模式将不同基因进行聚类，分析基因模块与特定性状之间的关系，筛选出有效生物标志物[5-6]。目前，WGCNA被广泛用于疾病、生理学、药物、进化和基因组的生物学研究[7]。

本研究拟通过对尿道下裂相关的基因表达数据分析整合后寻找出与尿道下裂发生相关的致病基因，为今后研究疾病的发病机制和靶向治疗提供线索。

1  材料与方法

1.1  数据获取

从GEO（https//www.ncbi.nlm.nih.gov/gds）数据库中下载与尿道下裂有关的芯片包含表达数据及临床信息（GSE35034）。其收录了来自日本国家儿童健康与发展研究所的尿道下裂儿童患者（n=23，中位年龄2.3岁）接受外科手术治疗的3种人类包皮成纤维细胞（human foreskin fibroblasts，hFFC）进行了全基因组筛选。用二甲基亚砜（dimethyl sulfoxide，DMSO）作为对照，分别暴露于10nmol/L双酚A（bisphenol A，BPA）、0.01nmol/L17β-雌二醇（estradiol，E2）和1nmol/L 2，3，7，8-四氯二苯并对二英（tetrachlorodibenzo-p-dioxin，TCDD）24h后測量来自尿道下裂患儿的3种hFFC的基因表达谱。所有样本均采用Agilent-028004 SurePrint G3 Human GE平台检测。

1.2  WGCNA共表达模块分析与功能富集分析

WGCNA是一种系统的生物学方法，常用于描述不同样本之间的基因关联模式[8]。使用“WGCNA” R包构建共表达网络[9]。选取合适软阈值，并筛选出权重最大模块。在WGCNA最佳模块中选择核心基因（R>0.8，P<0.01），通过利用R的“clusterprofiler”包对模块核心基因进行基因本体论（gene ontology，GO）和基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路功能富集，筛选出具有显著富集的功能和信号通路。

1.3  筛选差异基因与核心基因

利用R语言的“Limma”包对GSE35034所有转录数据进行差异分析，获取与尿道下裂发病相关的差异基因[10]。对E2、TCDD、BPA组数据筛选得到的差异表达基因取交集，获得表达一致的差异基因。将差异基因输入到String数据库中，筛选并绘制蛋白质互作网络（protein-protein interaction network， PPI）图[11]。使用Cytoscape软件对网络进行分析，选取连通性最高的前10个基因作为尿道下裂发病的枢纽基因。为使结果更为准确，将WGCNA核心模块基因、共同差异基因及PPI网络所筛选的基因，系统分析取其交集，筛选出重叠的关键基因。再对关键基因进行mRNA表达变化及受试者操作特征（receiver operating characteristic，ROC）曲线，计算ROC曲线下面积（area under curve，AUC），诊断验证。

1.4  统计学方法

本研究所有数据采用R软件分析，差异分析采用贝叶斯算法，相关性分析采用皮尔逊相关性分析，两组间的数据分析采用双t检验，适当的时候用t检验，以P<0.05为差异有统计学意义。

2  结果

2.1  样本预处理及基因共表达模块的构建

對原始数据进行预处理，共获得109 834个基因的表达数据。对样本进行离群性检验，然后进行聚类分析。首先对软阈值进行筛选，当软阈值β设置为8时，样本独立性较高，且可获得较好的连接关系。将基因动态切割，得到基因聚类树（图1）。再连锁分层聚类合并，最终得到15个模块，各模块基因个数从110～357不等。

2.2  基因模块与临床性状关联分析

计算模块的特征基因与样本来源的皮尔逊相关系数，并考虑到相关系数和P值，r=0.72，P=0.008作为具有临床意义的模块。然后对其进行基因重要性及基因表达水平分析，并绘制散点图分析模块内基因是否符合线性（图2）。

2.3  GO和KEGG富集分析

利用R的“clusterProfiler”包对选择模块基因进行GO和KEGG功能富集分析（图3）。GO富集结

果显示，模块基因主要参与肾脏发育，生殖系统发育，促进脂质生物合成等过程。在KEGG通路方面，主要富集于一些与细胞增殖分化等有关调节，如癌症通路，PI3K-Akt信号通路，蛋白质消化吸收，PPAR信号通路。

2.4  差异基因的PPI构建

从GEO数据库分别获得的E2、TCDD、BPA组数据，经过差异分析，得到满足条件的差异基因，再将三组数据取交集共得到共同差异表达基因93个，其中上调基因 37个，下调基因56个，并用火山图展示其表达的差异性（图4）。

将筛选出的共同差异基因导入STRING，筛选条件设置为：置信度≥0.15，互作最大值等于0。之后把STRING的计算结果导入Cytoscape 软件。应用cytoHubba插件筛选出节点打分排名前10的关键基因：CTNND2、SLITRK1、KCNF1、FAM107A、CTTNBP2、SCAMP5、TMEM132B、LAMP5、SYNDIG1、FBXL16。

2.5  关键基因表达验证和ROC曲线分析

为进一步验证核心分子准确性，对数据集中WGCNA核心模块基因、共同差异基因及Cytoscape软件所筛选的基因，整合分析，最后筛选出关键基因SYNDIG1和FBXL16为有效的生物标志物。

在数据集GSE121712中验证所筛选的基因在尿道下裂患者与正常组表达存在差异（图5），同时，ROC分析表明，发现SYNDIG1和FBXL16基因对于尿道下裂的存在具有较好的预测能力（图6）。

3  讨论

尿道下裂是一种男性常见生殖系统畸形，目前手术治疗效果不满意且术后并发症多，因此早期诊断预防具有重要意义。目前公认的尿道下裂发病可能与基因、环境等多因素相关，但其详细机制仍不明确。携带易感基因并存在环境内分泌干扰物暴露是其中一种较为合理的推论。本研究通过WGCNA方法，筛选与尿道下裂发病相关的基因，并进行了验证分析。

通过基因共表达网络，筛选出了与尿道下裂发病具有显著相关性的棕色模块。然后进行了功能富集分析。根据富集结果显示，关键基因在生物学过程方面主要参与肾脏发育，泌尿系统发育和泌尿生殖系统发育；而在功能和信号通路方面，主要参与一些与细胞增殖分化等有关调节，如： PI3K-Akt信号通路，卵巢类固醇生成，PPAR信号通路。其中，PI3K-Akt信号通路在先天性尿道下裂中发挥重要作用，PPAR信号通路与尿道发育相关[12-13]。相关研究表明，PI3K/Akt信号广泛参与调节各种细胞过程，包括增殖、分化、存活、迁移和凋亡等；PI3K/Akt的激活可以减少雄激素受体介导的转录[14]。

同时本研究进行差异分析和PPI网络构建，与棕色模块联合分析后，从中筛选出与尿道下裂生物发病过程最可能相关的两个基因SYNDIG1和FBXL16。SYNDIG1即突触分化诱导基因1，主要参与调控突触后受体密度[15]。

FBXL16即F-box和富含亮氨酸重复蛋白16 ，是一种F-box蛋白，参与几乎所有生命活动的调节，如细胞周期、增殖和凋亡。有研究表明，FBXL16可能与泛素化、磷酸化有关，可以促进细胞增殖[16]。Sato等[17]的实验也证明，FBXL16表达的缺失促进了细胞增殖。一项关于乳腺癌的研究中也发现FBXL16沉默抑制了细胞中p-AKT和p-mTOR的水平。FBXL16表达抑制乳腺癌上皮-间充质转化和血管生成也被证实。FBXL16还被证明是心血管祖细胞分化抑制因子之一[18-20]。在胚胎发育过程中泌尿系统及原始血管系统均起源于中胚层。原始血管的形成与血岛边缘细胞向血管内皮细胞的分化有关。男性尿道则形成主要依赖于尿生殖窦的中下段，生殖结节演变为阴茎，左、右尿生殖褶在中线闭合并形成尿道海绵体部。FBXL16已被证实可以抑制心血管祖细胞分化。而类似的机制也可能发生于尿道沟闭合、尿道形成时。其具体机制仍有待进一步探索和实验验证。

综上所述，本研究通过构建WGCNA构建共表达网络，分析处理后发现了SYNDIG1和FBXL16两个基因与环境共同作用后在尿道下裂生物发病过程中可能起着重要作用。其机制可能与PI3K/Akt/ mTOR信号通路相关。这些发现对于研究尿道发病机制，以及早期诊断、靶向治疗有重要参考价值。

利益冲突：所有作者均声明不存在利益冲突。

[参考文献]

[1] VAN D Z L F， VAN R I A， FEITZ W F， et al. Aetiology of hypospadias： A systematic review of genes and environment[J]. Hum Reprod Update， 2012， 18（3）： 260–283.

[2] MANSON J M， CARR M C. Molecular epidemiology of hypospadias： Review of genetic and environmental risk factors[J]. Birth Defects Res.， Part A， 2003， 67（10）： 825–836.

[3] WEI J， ZHANG J， WEI J， et al. Identification of AGXT2， SHMT1， and ACO2 as important biomarkers of acute kidney injury by WGCNA[J]. PLoS One， 2023， 18（2）： e0281439.

[4] ZHANG T， WONG G. Gene expression data analysis using hellinger correlation in weighted gene co-expression networks （WGCNA）[J]. Comput Struct Biotechnol J， 2022， 20： 3851–3863.

[5] MASON M J， FAN G， PLATH K， et al. Signed weighted gene co-expression network analysis of transcriptional regulation in murine embryonic stem cells[J]. BMC Genomics， 2009， 10： 327.

[6] LIANG W， SUN F， ZHAO Y， et al. Identification of susceptibility modules and genes for cardiovascular disease in diabetic patients using WGCNA analysis[J]. J Diabetes Res， 2020， 2020： 4178639.

[7] TIAN Z， HE W， TANG J， et al. Identification of important modules and biomarkers in breast cancer based on WGCNA[J]. Onco Targets Ther， 2020， 13： 6805–6817.

[8] WANG Y， ZHOU W， CHEN Y， et al. Identification of susceptibility modules and hub genes of osteoarthritis by WGCNA analysis[J]. Front Genet， 2022， 13： 1036156.

[9] LANGFELDER P， HORVATH S. WGCNA： An R package for weighted correlation network analysis[J]. BMC Bioinf， 2008， 9： 559.

[10] WEN J， REN T， ZHENG J， et al. Identification and verification of pivotal genes promoting the progression of atherosclerosis based on WGCNA[J]. Artif Cells， Nanomed， Biotechnol， 2023， 51（1）： 276–285.

[11] SZKLARCZYK D， GABLE A L， NASTOU K C， et al. The STRING database in 2021： customizable protein- protein networks， and functional characterization of user- uploaded gene/measurement sets[J]. Nucleic Acids Res， 2021， 49（D1）： D605–d612.

[12] CHAO C， JIAN Y Y， ZHAO X Y， et al. The involvement of hsa_circ_0000417 in the development of hypospadias by regulating AR[J]. Differentiation， 2020， 116： 9–15.

[13] MANSOUR M M， GOYAL H O， BRADEN T D， et al. Activation of penile proadipogenic peroxisome proliferator- activated receptor gamma with an estrogen： interaction with estrogen receptor alpha during postnatal development[J]. PPAR Res， 2008， 2008： 651419.

（下轉第18页）