生物信息学分析筛选骨关节炎关键基因的研究▲

蒋 捷 黄林科，4 韦林华，5 徐家科，3* 陈蔚蔚#

(1 广西医科大学第一附属医院创伤骨科手外科，广西南宁市 530021；2 广西医科大学再生医学研究中心，广西南宁市 530021；3 西澳大利亚大学生物医学与科学学院，澳大利亚珀斯 6009；4 广西医科大学第二附属医院骨科，广西南宁市 530006；5 广西医科大学附属南宁市传染病医院外科，广西南宁市 530023)

骨关节炎(osteoarthritis，OA)是影响中老年人生活质量最常见的慢性病之一，发病率为6%～32%[1]。OA是一种异质性疾病，具体病因不详，一般认为与遗传、年龄、肥胖及机械应力等有关[2]。OA病理变化主要包括4个方面：关节软骨退变、软骨下骨硬化、关节周围骨赘形成和滑膜炎症浸润[3]。对OA发病机制和治疗手段的研究，目前主要集中在软骨组织和软骨细胞上[4]，但研究结果并不理想。近年来，通过对滑膜组织及软骨下骨在OA发生发展过程中的作用的研究，得到了更多的认识[5]。滑膜炎引起的不正常滑膜液循环导致关节积液、引起关节疼痛，是早期OA的主要表现[6]。不同阶段OA伴有不同程度的滑膜炎症[7]。因此，滑膜炎症的研究对揭示OA的发病机制、疾病进展和治疗方案的选择有着重要的作用。

随着生物信息学的快速发展，越来越多的研究表明OA的发生发展与一系列的基因和信号通路相关[8]。利用GEO数据库(Gene Express Omnibus, http: //www.ncbi.nlm.nih.gov/geo)中关节滑膜基因表达谱芯片数据，筛选健康群体和OA患者滑膜组织的差异表达基因(differentially expressed genes，DEGs)，对DEGs进行聚类分析及功能富集分析，并进一步构建蛋白相互作用(protein-protein interaction, PPI)网络，筛选出与OA密切相关的关键基因，旨在进一步识别、鉴定OA发生发展中的靶基因，为OA发病分子机制的明确及分子标志物的筛选提供理论指导。现报告如下。

1 材料与方法

1.1 基因芯片数据 从GEO数据库下载基于GPL96[HG-U133A] Affymetrix Human Genome U133A Array平台的基因表达谱芯片数据集：GSE55457、GSE55235、GSE12021。3个数据集共包括了30例OA患者和29例健康对照者的滑膜组织(见表1)。所有数据集中的类风湿关节炎患者样本被排除，GSE12021数据集中基于GPL97平台的样本也被排除。

表1 3个GEO数据集信息

1.2 DEGs分析 利用GEO自带的在线工具GEO2R(https://www.ncbi.nlm. nih.gov/geo/geo2r/)对每个数据集的OA患者及健康对照者的滑膜组织进行DEGs分析，按照调整后P值(adjustedP)<0. 05，|log2FC|>1(FC表示DEGs上调或下调倍数)的筛选条件分别筛选3个数据集中OA和健康对照者滑膜组织的DEGs。使用在线分析工具jvenn(http://jvenn.toulouse.inra.fr/app/index.html)对3个数据集中共同的DEGs进行筛选，确定DEGs。

1.3 DEGs的GO功能富集及KEGG通路富集分析 利用在线分析工具DAVID 6.8(https://david.ncifcrf.gov/)[9]对筛选的DEGs进行GO(Gene Ontology, GO)功能富集分析和KEGG(Kyoto Encyclopedia of Genes and Genomes,KEGG)通路富集分析，以P<0.05为差异有统计学意义。

1.4 差异基因PPI网络的构建与关键基因分析 应用在线分析数据库STRING(http://string-db.org/cgi/input.pl)[10]，将筛选出的DEGs导入数据库，进行差异基因PPI分析。利用Cytoscape软件构建PPI网络[11]，并应用其中的Cytohubba插件选取前10 位的差异基因作为蛋白质网络中的hub基因。MCODE插件对PPI进行模块化分析(参数：degree cutoff=2, node score cutoff=0.2, K-core=2, max. cepth=100)，选取score>3且基因数大于4个的子网络作为显著性模块。

2 结 果

2.1 健康对照者与OA患者滑膜组织的DEGs 分别从基因芯片数据集GSE55457、GSE55235、GSE12021中筛选出786、1 048、808个差异基因。在3个数据集中均有的DEGs为126个(图1)。

图1 GSE55457、GSE55235、GSE12021的DEGs

2.2 DEGs的GO功能富集和KEGG通路富集分析结果 GO富集分析表明，滑膜组织中DEGs参与的主要生物学过程有：细胞对促肾上腺皮质激素释放激素刺激的反应、通过mRNA剪接体对mRNA的剪接、负向调控基因表达、正向调控单核细胞聚集、正向调控成纤维细胞增殖、蛋白转运、调控炎症反应、骨骼肌细胞分化、调节呼吸过程中气体交换、正向调控细胞迁移；细胞成分包括核质、细胞核、生长锥膜；分子功能主要集中于核苷酸结合、肝素结合、RNA聚合酶Ⅱ转录因子活性与配体激活序列特异性DNA结合、类固醇激素受体活性。详见表2。

表2 DEGs的GO分析结果

KEGG富集通路分析显示，滑膜组织DEGs的主要信号通路富集在包括MAPK信号通路、EB(Epstein-Barr，EB)病毒分子感染通路、Ⅰ型人类T细胞白血病病毒分子感染通路、癌症中蛋白聚糖分子通路、恶性肿瘤分子通路、破骨细胞分化通路、PI3K/Akt信号通路、胰岛素信号转导通路、肿瘤中miRNAs分子通路、ErbB信号通路、HIF-1信号通路和膀胱癌分子相关通路。详见图2。

图2 DEGs参与的KEGG信号通路

2.3 DEGs的PPI网络分析 将DEGs导入STRING和Cytoscape构建出PPI网络，该网络由84个结点和233条边构成(图3)。Cytohubba插件确认的度值前10的hub基因是(图4)：JUN、BTG2、ATF3、DUSP1、HNRNPA1、SF1、SRSF7、TRA2B、TRA2A、PPP1R15A。

图3 DEGs编码PPI网络

图4 hub基因及其扩展PPI网络

MCODE插件分析出score>3且基因数多于4个的子网络有4个(图5)，其中最显著的相互作用模块(score:7)由7个结点和21条边构成，包含了TRA2B、SRRM2、SRSF11、TRA2A、HNRNPA1、SRSF7、SF1 7个hub基因，这些基因存在非常复杂的相互作用。

图5 MCODE确认的网络核心模块

3 讨 论

OA是世界范围内致残的主要因素[12]。但OA发病的病理机制仍然不明确，因此了解OA的发病机制将有助于早期诊断、判断预后和发现药物的新靶点[13]。

本研究从GEO数据库获取健康对照者和OA患者滑膜组织的基因表达谱芯片数据，利用一系列在线分析工具(GEO2R、jvenn、DAVID、STRING、Cytoscape)，结合生物信息学方法对二者的DEGs进行分析，筛选出了126个DEGs。对DEGs进行GO功能富集，这些基因主要参与细胞对促肾上腺皮质激素释放激素刺激的反应、通过mRNA剪接体对mRNA的剪接、负向调控基因表达、正向调控单核细胞聚集、正向调控成纤维细胞增殖等生物学过程。通过KEGG通路富集分析，我们确定了与本团队前期研究相一致的影响OA发病机制通路，包括MAPK、PI3K/Akt和HIF-1信号通路。研究表明，抑制MAPK信号通路可以减轻滑膜炎症，改善炎症诱导的软骨细胞凋亡，促进软骨再生[14-15]。此外，PI3K/Akt信号通路对细胞的增殖、凋亡过程具有重要的调节作用，在OA的发生、发展过程中起着一定的作用。有研究报道PI3K/Akt信号通路在OA中高表达[16]，抑制PI3K/Akt信号通路会降低胶原诱导关节炎小鼠滑膜钙粘蛋白-11的表达和大鼠关节炎的炎症反应[17-18]。HIF-1是细胞对缺氧反应的主要转录调节因子，由于软骨组织生长在一个低氧的微环境中，HIF-1在软骨稳态的平衡维持中发挥重要作用，活化的HIF-1与软骨下骨重构、新生血管生成有密切关系[19-20]。

同样，我们构建了DEGs的PPI网络，利用Cytoscape软件中的Cytohubba插件对差异基因PPI网络中连接最为紧密且位于重要结点位置的基因进行分组，筛选出了JUN、BTG2、ATF3、DUSP1、HNRNPA1、SF1、SRSF7、TRA2B、TRA2A、PPP1R15A等hub基因。JUN是编码激活蛋白-1的亚单位，白细胞介素-1能够通过激活JUN诱导基质金属蛋白酶-13产生，从而导致软骨基质降解[21-22]。JUN编码的c-JUN蛋白对调节软骨细胞凋亡起着重要的作用[23]。ATF3在OA软骨中的表达显著升高，ATF3可以调节炎症因子如白介素细胞-6的分泌，影响软骨的代谢[24]。TRA2B是一种重要的RNA结合蛋白，在预后较差的骨肉瘤患者中高表达，miR-206可下调TRA2B表达诱导骨肉瘤细胞凋亡[25]。BTG2属于BTG/TOB家族成员，可有效下调PI3K和Akt磷酸化抑制骨肉瘤细胞迁移侵袭能力[26]。敲除BTG2导致小鼠的腰椎后移，与BMP信号通路活化减弱相关[27-28]。HNRNPA1是一种RNA结合蛋白，参与信使RNA的代谢和转运，突变可能引起多系统蛋白病。TRAF6调节HNRNPA1的泛素化失衡会诱发骨髓异常增生综合征[29-30]。SRSF7、SRRM2、SRSF11、BTG2、SF1、HNRNPA1、PPP1R15A暂时未见与OA相关的报道，可能是新的潜在的OA发生相关基因。

在本研究中，我们同时发现了一些与OA有关的新的标志物。DUSP1在3个数据集的OA患者中相对健康对照者表达量明显减少。DUSP1是丝裂原活化蛋白激酶(mitogen-activated protein kinases, MAPKs)的内源性抑制物，下调或者敲除DUSP1可以增强MAPKs作用[31-32]。在OA中，白细胞介素-1β诱导的MAPKs的活化导致炎症因子诱导的软骨细胞凋亡，基质金属蛋白酶-3、基质金属蛋白酶-13等酶的大量产生导致软骨基质降解[33-34]。因此，是否可以通过增加DUSP1的表达，减少MAPKs的活化，从而达到减轻炎症、保护软骨、延缓OA发生的作用，值得进一步研究。

本研究的不足之处：(1)未能进一步做到hub基因的GO和KEGG分析；(2)未在实例标本中证明hub基因的表达情况；(3)未能分析MCODE插件中得出的各子网络的作用。在今后的研究中，我们会完善相关实验，得出更完整的结论。

综上所述，在对本研究筛选得到的DEGs如JUN、ATF3、DUSP1等这些分子生物功能及信号通路的进一步研究中，对OA的诊治将会有更多的新发现。BTG2、HNRNPA1、SF1、SRSF7、TRA2B、TRA2A、PPP1R15A这些基因可能是潜在的影响OA的发生发展的新基因，虽然其在OA中的价值仍然需要通过大量的实验证实，但为今后的进一步研究提供了重要的线索和参考信息。