原发性骨质疏松症相关核心基因的生物信息学分析

杨洋 薛建华 虞俊波 顾琪琪 张亚峰 刘佳佳*

1. 南通大学附属医院创伤中心，江苏 南通 226001

2. 南通大学附属医院骨科，江苏 南通 226001

原发性骨质疏松症(osteoporosis，OP)是一种以骨量减少、骨的微结构退化，骨骼脆性增加和骨折发生风险增加为特征的全身性骨代谢疾病[1]。OP常无明显的临床症状，一般在骨折发生后才被发现，主要累及脊柱、髋部、腕部和肱骨近端等部位，故OP又被称为“无声的浩劫”。因此，及时地预测OP对防治骨质疏松骨折尤为重要[2]。

骨代谢的平衡由成骨细胞的骨形成作用和破骨细胞的骨吸收作用共同维持，两者之间的动态平衡共同维持骨骼微环境的稳态[3]。既往研究表明，过度激活的破骨细胞以及功能降低的成骨细胞是OP的重要病因。针对Wnt、NF-κB、MAPK等通路设计开发药物，重塑骨代谢平衡，在OP的治疗中具有重要意义[4-5]。然而，现有抗OP药物如双膦酸盐、特立帕肽等，仍具有非典型性骨折、下颌骨坏死、费用昂贵等特点[6]。探究OP发病机制，寻找新的药物靶点，在OP的临床治疗中具有重要意义[7-8]。

因此，为进一步探索OP的分子机制，本研究拟通过分析基因表达数据库(gene expression omnibus，GEO)提供的OP相关基因芯片数据，筛选差异表达基因，并对这些基因进行核心基因和信号通路的生物信息学分析，为进一步在分子水平研究OP的发病机制提供理论依据。

1 材料和方法

1.1 一般资料

研究号 GSE35958 基因表达谱芯片数据来源于美国国立生物技术信息中心的公共基因表达数据库(gene expression omnibus ，GEO)(https://www.ncbi.nlm.nih.gov/geo/)。该数据集由Benisch 等[9]于2012 年发表，由HGU133 Plus 2.0 ( Affymetrix公司,美国)注释，芯片平台为 GPL570。该数据集包含2组来自不同人群的人股骨头骨髓间充质干细胞(HMSC)，共计9例样本。其中4例样本的供体为行全髋关节置换术的非OP老年患者，另5例样本的供体为股骨颈低能量骨折伴有OP的老年患者。

1.2 差异基因的数据分析

利用R软件对芯片原始数据进行注释和过滤，使用Bioconductor (http://www.bioconductor.org)提供的RMA算法对各原始芯片数据进行背景校正及归一化等预处理，得到标准化的数据。选取校正后的P<0.05且∣logFC∣≥2为阈值，获得该数据集的差异表达基因(differentially expressed genes，DEGs)。

1.3 GO功能和KEGG信号通路富集分析

本研究采用Metascape(http://metascape.org/gp/)在线分析工具进行差异基因富集分析。Metascape整合多个数据库资源(GO、KEGG、UniProt和DrugBank等)，本研究使用的是基因本体( gene ontology，GO)及京都基因和基因组百科全书( kyoto encyclopedia of genes and genomes，KEGG) 两个数据库，分别对差异基因进行GO功能富集以分析其参与的生物学进程和KEGG 信号通路富集分析。

1.4 PPI与核心模块分析

通过STRING(search tool for the retrieval of interacting genes /proteins) 11.0在线工具对差异基因进行蛋白质-蛋白质相互作用( protein-protein interaction，PPI) 网络分析。将STRING的计算结果导入Cytoscape 3.6.1软件，利用插件MCODE和Cytohubba分别寻找蛋白相互作用网络中与OP相关的核心基因和核心模块。

2 结果

2.1 差异基因的获取

本研究选用基因表达谱芯片GSE35958，通过R软件采用RMA算法对数据进行预处理后，根据筛选标准筛选出569个差异基因，其中上调7个，下调562个。通过ImageGP(http://www.ehbio.com/ImageGP/)在线工具绘制OP样本与非OP样本差异表达基因热图。热图的颜色反映基因表达量的水平，白色代表升高，蓝色代表降低，颜色的亮度代表基因升高及降低的程度，见图 1。

2.2 基因富集分析

将筛选出的差异基因导入Metascape后，分别勾选GO中的生物学过程(biological processes)、细胞定位(cellular components)、分子功能(molecular functions)和KEGG中的信号通路(pathway)提交分析。可以看出，差异基因功能主要富集在mRNA前体内含子结合、核体、组蛋白修饰、mRNA 3’端加工和调控结合等(图2)，通路主要富集在叶酸生物合成、粘合连接和致病性大肠杆菌感染等(图3)。然而，这3条信号通路中并不包含本研究筛选出的核心基因。基于此，本研究继续对10个核心基因再次进行KEEG通路富集分析，得到一条核心基因显著富集的信号通路，即泛素介导蛋白水解信号通路(hsa04120)(图4)。

2.3 PPI构建及分析

利用STRING在线工具对569个差异表达基因进行 PPI构建，设定最低相互作用分数为 0.8，得到了由 517个节点、363条连线构成的相互作用关系图(图5)。将STRING 的计算结果导入Cytoscape 3.6.1软件后，采用MCC算法得到前10个核心基因，按值的大小由高到低分别是TCEB1、CUL2、KBTBD6、KBTBD7、ASB8、KLHL42、ASB5、FBXO11、ANAPC10、CDC23，且这10个基因全部为下调的差异表达基因。此外，用MCODE插件对该蛋白相互作用网络进行分析，得到了由上述10个核心基因构成的核心模块(图6)。

3 讨论

骨质疏松症是一种年龄相关的慢性骨骼代谢疾病，主要以骨量减少、骨微结构改变和骨脆性增加为特点。目前，OP的发病机制尚未完全阐明，但目前已知OP的发生具有一定的遗传因素，与特定基因的表达有关[10-11]。

本研究分析了2组来自不同人群的人股骨头骨髓间充质干细胞(HMSC)，共计9例样本，通过GEO数据库中的基因表达谱芯片GSE35958获得差异基因，并分析差异基因GO功能富集和KEGG通路富集，发现泛素介导蛋白水解信号通路显著激活，为OP发病机制的研究提供了理论依据和研究方向。在本研究中，KEGG 最终筛选出的富集通路为泛素介导蛋白水解信号通路，该通路包含CUL2、CDC23、ANAPC10 3个核心基因。而在骨代谢过程中，破骨细胞首先粘附在骨组织，然后通过酸化和蛋白水解完成骨吸收[12]，提示该通路可能与OP有关联。另外，就分析得到的核心基因而言，现有研究主要集中在肿瘤领域，而本研究为这些基因参与OP的发生提供了一定依据，有助于拓展核心基因的功能。

CUL2编码的蛋白质为Cullin蛋白家族成员，包括Elongin B和C，Rbx1和各种底物识别受体的支架蛋白，形成E3泛素连接酶[13-14]。有学者发现，E3泛素连接酶可以与HIFβ结合，促进HIFα的泛素化和蛋白酶降解[15-16]，而HIFα又可与HIF1β结合形成HIF转录因子，调节基因表达[17]。HIF下游靶基因中包括血管内皮生长因子A(VEGFA)，可促进细胞生长和微血管的增殖[15]。以上结果提示CUL2编码的蛋白质可能会通过调控HIF表达以影响骨血管再生和成骨细胞活性，从而影响OP的发生发展。

TCEB1编码的蛋白质Elongin C为转录因子B(SIII)复合物的一个亚基，SIII复合物由Elongin A / A2、B和C组成。其中，Elongin B和C是调节亚基，其支架蛋白为CUL2编码的蛋白质。因此，TCEB1可能与上述的HIF调控下游靶基因VEGFA有关，具体结果有待进一步探讨。

CDC23编码的蛋白质与酿酒酵母Cdc23有很强的相似性，后者为细胞周期G2/M转换所必需的蛋白质。该蛋白是后期促进因子(APC)的组成部分，APC由8个蛋白亚基组成，在真核细胞中高度保守，可催化细胞周期蛋白的形成。有研究发现，CDC23在甲状腺癌组织中过度表达，而在正常和增生的甲状腺组织中缺失，并且CDC23为细胞周期和细胞生长的关键调节因子[18]。因此，验证 CDC23是否参与OP的骨偶联调控值得深入研究。

ANAPC10编码的蛋白质为后期促进因子(APC)的一个亚基，由185个氨基酸组成，具有一个保守的核心，是细胞周期调控的泛素蛋白连接酶，在细胞从有丝分裂到G1的进程中起着至关重要的作用[19-20]。由于APC10与SP1蛋白存在功能相关，而SP1蛋白在乳腺癌、甲状腺癌和胰腺癌等癌组织中与正常组织相比呈现高表达[21-22]。因此，APC10的失调可能与肿瘤的发生有关。但是，ANAPC10与OP是否存在相关性，目前尚无证据支撑，可作为研究的一个方向。

综上所述，本研究分析得到的10个核心基因和一条信号通路为研究OP 的发生机制和药物治疗提供了新的思路，或为OP药物研发提供作用靶点。然而，本研究所得到的结论尚需通过基础实验来作进一步验证。