基于GEO 数据库分析膝骨关节炎软骨的关键差异基因

许展仪，朱根福

（广州中医药大学第三附属医院，广东 广州）

0 引言

骨性关节炎是一种以软骨退变、滑膜炎症、骨赘形成和软骨下骨硬化为特征的慢性关节疾病。其典型的体征和症状包括疼痛、肿胀和僵硬，常伴有功能减退和活动受限[1]。这是一种进展缓慢的致残性关节紊乱，显著降低了生活质量。我国60岁及以上人群膝骨性关节炎的患病率为35%-50%[2]。尽管对OA的形成和发展的机制和病因进行了广泛的研究，但OA的病因仍不清楚。流行病学研究表明，OA是一种复杂的多基因疾病，有许多环境和遗传危险因素，其中导致疾病进展的因素之一是遗传成分[3]。在过去的15年里，研究的重点是寻找骨关节炎的易感部位。基因芯片筛选出的差异表达基因，可用作早期诊断和靶向治疗的生物标记物[4]。目前，随着高通量测序技术的发展，人们对骨关节炎基因表达谱进行了大量研究，筛选出数千个差异表达基因。然而，由于样品的异质性或测序平台的不同，表达的mRNA的结果与不同的基因谱不一致。因此，在OA中还没有确定可靠的结果。然而，集成的生物信息学方法将解决这些缺点，并确定参与OA的hub基因。

在本研究中，我们下载了微阵列数据集GSE117999[5]，并筛选出膝关节OA患者和正常对照组软骨之间的差异表达基因（DEGs）。应用DEGS的GO和KEGG富集分析，构建了蛋白质-蛋白质相互作用(PPI)网络的功能模块分析。本研究旨在鉴定hub基因，探索OA发生的内在分子机制。

1 材料和方法

1.1 基因芯片原始数据骨性关节炎软骨组织基因芯片数据集为GSE117999，通过GEO数据库（https://www.ncbi.nlm.nih.gov/GEO/）[6]下载得到，该数据集的平台为：Agilent-072363 SurePrint G3 Human GE v3 8x60K Microarray 039494 [Feature Number Version]，该数据集为转录组表达芯片，共包括12例骨关节炎患者（GSM3316970、GSM3316972无表达数据）和12例无骨关节炎行关节镜部分半月板切除术患者（GSM3316976、GSM3316980无表达数据）的软骨组织。

1.2 原始数据预处理和差异表达基因（DEGs）筛选 利用GEO2R在线分析工具[7]对数据集GSE117999进行分析，其借助R语言（R3.2.3）软件中Limma 3.2.8工具包对数据进行差异表达分析，筛选出差异表达基因，筛选条件设置为P<0.05，差异倍数（Fold Change，FC）≥1。利用DAVID中“Gene ID Conversion”对基因探针与基因库中的基因名进行匹配。

1.3 差异表达基因GO和KEGG富集分析 借助DAVID在线数据库（https://david.ncifcrf.gov/）进行GO富集分析[8]，其包括生物过程（BP）、细胞成分（CC）和分子功能（MF）3个方面的功能分析，以了解其相关功能。与此同时，本研究还利用KEGG在线数据库（https://www.genome.jp/kegg/pathway.html）对差异表达基因进行KEGG信号通路分析[9]，以了解骨关节炎软骨病变中潜在的KEGG信号通路，P<0.05被认为差异具有统计学意义。

1.4 蛋白互作网络（PPI）构建和hub 基因筛选

通过STRING数据库（https://string-db.org/）分析蛋白互作关系，该在线数据库包含了5090个物种的2458万蛋白质[10]。借助可视化软件Cytoscape 3.7.0（http://www.cytoscape.org/）构建蛋白互作网络[11]，使用CytoHubba插件进行PPI网络的hub基因分析，通过MCC、MNC、Degree这3种不同的算法构建3个蛋白互作网络，得分最高的前10个节点作为hub基因。

2 结果

2.1 差异表达基因筛选结果

与正常组相比，在膝骨关节炎软骨中共有85个DEGs，其中9个未能找到对应的基因名称，最终得到76个差异表达基因。其中上调的DEGs有33个，差异上调最明显的包括ATP结合盒式亚家族B成员5（ATP Binding Cassette Subfamily B Member 5，ABCB5）、蛋白磷酸酶 6调节亚基 2（Protein Phosphatase 6 Regulatory Subunit 2，PPP6R2）、防 御 素 α3（Defensin Alpha 3，DEFA3）、精子尾部外密纤维 2样（Outer Dense Fiber Of Sperm Tails 2 Like，ODF2L）、溶质载体家族 3成员 1（Solute Carrier Family 3 Member 1，SLC3A1）；下调的 DEGs有 43个，差异下调最明显的包括重组人序列相似性家族224成员A（Family With Sequence Similarity 224 Member A，FAM224A）、母系表达3（Maternally Expressed 3，MEG3）、软骨素样蛋白 2（Chordin Like 2，CHRDL2）、基质金属肽酶 3（Matrix Metallopeptidase 3，MMP3）等。见图1。

2.2 基因功能富集分析

差异表达基因GO富集结果显示，在分子功能中，差异表达基因主要涉及MHC II类受体活性、ATP酶活性；在细胞成分中，差异表达基因主要集中于反高尔基网络膜、MHC II类蛋白复合物、内质网膜腔侧的组成部分、核体；在生物过程中，差异表达基因主要与防御真菌、抗菌体液反应、对革兰氏阳性细菌的防御反应、慢性炎症反应、先天免疫反应有关。见表1。KEGG信号通路富集结果显示，差异表达基因主要富集于产生IgA的肠道免疫网络、类风湿关节炎、哮喘、结核、移植物抗宿主病、同种异体移植排斥等信号通路。见表1。

表1 骨关节炎患者中差异表达基因本体、信号通路富集分析结果

2.3 PPI 蛋白互作网络和hub 基因分析

分析结果通过3种不同的算法构建了3个PPI网络，得到了相似的结果。Cathelicidin抗菌肽（Cathelicidin Antimicrobial Peptide，CAMP）、防御素 α3（Defensin Alpha 3，DEFA3）、防御素α4（Defensin Alpha 4，DEFA4）在PPI网络中的各项算法评分是最高的，并且这3个hub基因的表达都是上调的。见图2。

图2 骨关节炎软骨中差异表达基因的蛋白质-蛋白质网络图

3 讨论

骨关节炎是世界上最常见的退行性疾病，影响大小关节。骨关节炎通常影响整个关节组织，包括滑膜、软骨下骨和软骨[12]。滑膜炎是骨关节炎的共同特征，即使在早期疾病中也是如此。滑膜增生和组织肥大在晚期骨关节炎中有重要意义。越来越多的证据表明滑膜炎在骨性关节炎的发生和发展中起着关键作用。进一步了解骨关节炎相关滑膜炎的分子和细胞变异性可以为关节炎的病因提供研究方向[13]。GEO数据库基因芯片研究有助于确定遗传因素在骨关节炎风险中的重要作用。了解影响骨关节炎发生和发展的关键基因将对于疾病的早期治疗的发展是必要的。

本研究中，基因表达谱分析揭示了与骨关节炎相关的hub基因和通路，并使我们能够确定治疗策略的靶点。应用生物信息学方法对原始数据进行分析，与正常组相比，我们在骨关节炎软骨中识别出76个DEG，其中上调的DEGS有33个，下调的DEGS有43个。利用DAVID在线数据库分析差异表达基因的GO和KEGG富集结果，其主要涉及MHC II类受体活性、ATP酶活性；主要集中于反高尔基网络膜、MHC II类蛋白复合物、内质网膜腔侧的组成部分、核体；并且显著富集于产生IgA的肠道免疫网络、类风湿关节炎、哮喘、结核、移植物抗宿主病、同种异体移植排斥等信号通路。在关节炎组织中还激活了免疫途径，即同种异体移植排斥以及产生IgA的肠道免疫网络，表明免疫反应参与了这些疾病的发生。炎症和细胞因子在RA和某些OA病例中起重要作用[14]。先天免疫反应与软骨的损伤存在着一定的关系，彼此的关系对OA的发展产生影响[15]。TLR增加IL-1β 和基质金属蛋白酶进而诱导软骨细胞的炎症反应，加重OA的恶化程度[16]。

关节软骨是一种结缔组织，主要由细胞外基质和透明软骨细胞组成，细胞外基质包括蛋白聚糖和胶原蛋白。正常软骨细胞外基质的合成和降解处于一种动态平衡的状态，如果基因表达受到干扰并且平衡被打破时，将发生OA。本研究借助STRING在线数据库和Cytoscape软件分析获得了前10个HUB基因：CAMP、DEFA3、DEFA4、S100A8、ARGLU1、MMP3、TOLLIP、MSH2、ESR1、DUT。研究这些hub基因可能有助于膝骨关节炎的早期诊断和治疗。基质金属蛋白酶（MMPS）是一组内肽酶，可以降解细胞外基质蛋白并具有同源结构，参与软骨破坏过程，从而参与骨关节炎的发生。MMP-3作为基质金属蛋白酶家族的重要成员之一，主要由软骨细胞、滑膜细胞分泌，不仅参与细胞外除多糖基质以外的所有基质的降解，还参与激活Ⅱ型胶原间质胶原酶的降解，从而推动软骨的破坏，促进关节炎的进程[17-18]。冯震[19]通过收集患者血液样本分析得知雌激素受体α(ESR1)内的rs2228480和rs2234693基因多态性与膝OA风险增加有关；rs9340799和rs2228480之间的相互作用以及含有rs2228480-A和rs2234693-C等位基因的单倍型与膝OA风险增加相关。Toll作用蛋白（Toll-interacting protein, Tollip）是一种接头蛋白，可在Toll样受体（Toll like receptors, TLRs）信号通路中起负反馈调节作用，TLR被广泛的病原体、细胞因子或其他特定分子激活，从而触发先天免疫反应，进而参与骨关节炎的进程[20]。IL-1刺激增加了软骨细胞S100a8和S100a9 mRNA和蛋白质水平，S100a8和S100a9的软骨细胞mRNA表达在早期小鼠OA 中升高，而在晚期OA中则没有。然而，与小鼠炎症性关节炎中的软骨细胞中的S100A8和S100A9两者的阳性反应相反，随着小鼠OA的进展，软骨细胞中的S100A8染色消失。S100A8和S100A9在小鼠和人类骨关节炎期间积极参与滑膜激活和关节破坏的调节[21]。Neidlin M等通过加权基因共表达网络分析（WGCNA）鉴定在四个关节组织中保存的基因模块发现8个基因（CSN1S1、APOD、FAP、COL5A2、MXRA5、DEFA3、DEFA4、S100A8）在其中 3个关节组织中差异表达[22]。本研究通过cytohubba插件分析发现一些新的与OA有关的基因，如DEFA3、DEFA4、MSH2等，可能作为 OA 潜在诊断和治疗的靶点。

通过生物信息学对GEO数据库关于骨关节炎数据集的分析，有助于确定骨关节炎发病机制中的有用靶点和途径。CAMP、DEFA3、DEFA4、S100A8、ARGLU1和MMP3等基因可能是调控骨关节炎发生发展的重要基因，这可作为骨关节炎诊断和治疗靶点的研究提供新的思路的方向。本研究的局限性在于仅在生物信息学方面筛选的差异表达基因和分析分子机制，其需要进一步实验研究加以证实。