摘要
骨关节炎（Osteoarthritis，OA）导致的疼痛和残疾已经极大地影响了患者的生活质量。近期许多研究表明，骨关节炎与糖酵解之间有密切关联。然而，糖酵解对OA的影响在很大程度上仍不清楚。为此，通过生物信息学分析鉴定和验证可能参与OA滑膜炎的糖酵解相关基因。从Gene Expression Omnibus（GEO）数据库中选择微阵列数据集GSE55235。通过维恩图和NetworkAnalyst筛选糖酵解表型差异基因。通过GO和KEGG富集分析以及蛋白质-蛋白质相互作用（PPI）分析筛选出表型差异基因。此外，利用Cytoscape软件和数据库STRING构建枢纽基因网络。然后通过分析GSE206848数据库上的受试者工作特征（ROC）曲线，确认了枢纽基因。考虑到糖酵解与免疫的关系，应用CIBERSORTx分析了OA中的免疫浸润。研究发现了4个基因，包括表皮生长因子受体（Epidermal growth factor receptor，EGFR）、血管内皮生长因子A（Vascular endothelial growth factor A，VEGFA）、Egl-9家族缺氧诱导因子3（Egl nine homolog 3，EGLN3）和DNA损伤诱导转录子4（DNA Damage Inducible Transcript 4，DDIT4），为枢纽基因。ROC分析表明，几乎所有枢纽基因在GSE206848中都具有良好的诊断特性。本研究发现4个表型差异基因是OA期间滑膜炎潜在的诊断生物标志物和治疗靶点，这也许能完善OA的发病机制。
Abstract
Osteoarthritis significantly compromises patients' quality of life due to the associated pain and disability. A growing body of evidence underscores a robust association between osteoarthritis (OA) and glycolysis, yet its influence on OA largely remains to be elucidated. This research aimed to identify and authenticate glycolysis-associated genes potentially implicated in OA synovitis through a comprehensive bioinformatic approach. The microarray dataset GSE55235 was extracted from the Gene Expression Omnibus(GEO) database. Glycolytic phenotype-specific differentially expressed genes were identified via a Venn diagram and NetworkAnalyst. Phenotypic variants were further screened using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, in addition to protein-protein interaction (PPI) analysis. A hub gene network was established utilizing Cytoscape software and the STRING database. Hub genes were pinpointed via Receiver Operating Characteristic (ROC) curve analysis of the GSE206848 database. Given the established interplay between glycolysis and immunity, this investigation employed CIBERSORTx to assess immune infiltration in OA. This study discerned four genes, namely EGFR, VEGFA, EGLN3, and DDIT4, as critical to the investigation. ROC analysis confirmed that the majority of these key genes exhibited reliable diagnostic attributes in the GSE206848 dataset. This research unveiled four phenotypically unique genes as potential diagnostic markers and therapeutic targets for synovitis in OA, thereby offering new insights into the disease's pathogenesis.
Keywords
骨关节炎是一种对老年人影响最大的退行性关节疾病[1]。它也是导致残疾的主要原因，其特点是软骨退化、滑膜炎症和软骨下硬化伴骨赘形成等关节结构的病理变化[2-3]。骨关节炎的危险因素有很多，例如陈旧性关节损伤，肢体或关节发育不良与异常，从事损伤关节的相关工作或者遗传等[4]。现在随着肥胖与人口老龄化问题逐渐显现，OA的发病率也在不断升高。然而，骨关节炎并不表现出早期症状，现有的方法无法检测出其微小的变化，所以需要依靠灵敏度更高的生物标志物来帮助OA的早期诊断，以此来促进完善OA的发病机制研究。
糖酵解是指将葡萄糖或糖原分解为丙酮酸，ATP和NADH+H的过程，此过程中伴有少量ATP的生成。这一过程是在细胞质中进行，不需要氧气，每一反应步骤基本都由特异的酶催化。最近，糖酵解也成为了各种疾病研究的热点，例如癌症、糖尿病、遗传病等。
长久以来，研究者们一直认为OA主要是由于机械因素[5]。然而，越来越多的证据表明，代谢和炎症过程也是引起OA的关键原因。有研究发现，在炎症环境中，糖酵解的途径会增强，糖酵解代谢紊乱可导致软骨细胞肥大和细胞外基质降解，进而加快OA的进展[6]。由此可见，调节糖酵解代谢可能对于预防和治疗OA具有重要意义[7]。然而现在大多数研究都集中于糖酵解与关节软骨的关系，目前对于骨关节炎与滑膜的研究却很少，为了完善OA的发病机制理论，糖酵解代谢与滑膜组织的联系也需要更深入地研究与探索。本研究从基因表达综合数据库（GEO）中获取OA病例和健康受试者的数据，用于生物信息学分析，目的是筛选差异表达基因（DEG）。然后，将识别出的DEG与糖酵解数据集相交，以获得关键的GDEG。通过构建蛋白质-蛋白质相互作用（PPI）网络，本研究检测到某些中心GDEGs是OA发生的关键调控因子，是潜在的靶点。然后，对这些关键的GDEG进行GO与KEGG分析。此后，分析GSE206848数据并绘制ROC曲线，验证关键基因是否具有良好的诊断特性。最后，考虑到糖酵解与免疫力之间的潜在关系，使用CIBERSORTx分析了OA中的免疫浸润。总之建立一种与糖酵解相关的特异基因可为探讨OA机制以及OA的诊断和治疗提供更多见解。
1 材料与方法
1.1 收集材料
从GEO数据库（https://www.ncbi.nim.nih.gov/geo）中下载基于GPL96[HG-U133A] Affymetrix Human Genome U133A Array检测平台的基因微阵列数据集GSE55235，该数据集共包含30例滑膜样本，选择其中10例健康人为对照组，10例骨关节炎为实验组。
1.2 糖酵解相关差异基因筛选
使用NetworkAnalyst（https://www.networkanalyst.ca）对数据集进行分析，筛选条件设定为P＜0.05，|Log2FC|≥1，筛选OA滑膜组织和正常滑膜组织之间的差异表达基因（DEGs），并删除重复项。糖酵解相关基因从GSEA数据库（https://www.gsea-msigdb.org）获取，通过在线韦恩图获取糖酵解相关差异表达基因，通过Hiplot Pro在线平台（https://hiplot.com.cn）绘制糖酵解相关差异表达基因表达热图和相关性热图。
1.3 差异表达基因的GO和KEGG富集分析
利用SangerBox在线数据库（http://sangerbox.com）对筛选出的糖酵解相关差异基因进行基因本体功能富集分析（GO）和京都基因与基因组百科全书通路富集分析（KEGG），GO分析包括生物学过程、细胞成分和分子功能。FDR方法选择Benjamini & Hochberg，设定P＜0.05，FDR＜0.1，对结果进行可视化。
1.4 PPI网络构建与核心基因筛选
将上述糖酵解相关差异基因上传到STRING数据库（https://string-db.org），设置物种为智人，蛋白最低相互作用分值为0.4。将节点信息以TSV格式保存并导入Cytoscape3.9.1，构建蛋白-蛋白相互作用网络（PPI），利用cytoHubba插件按照degree值筛选前10个核心（Hub）基因。通过Cytoscape对Hub基因参与的部分信号通路进行结果可视化。
1.5 TF-miRNA-Hub基因调控网络构建
通过miRWalk（Version 3.0，http://mirwalk.umm.uni-heidelberg.de/）预测Hub基因相应的miRNA，并在Targetscan和miRDBand数据库进行验证。通过TRRUST（Version 3.0，https://www.grnpedia.org/trrust）数据库预测Hub基因相应的转录因子（Transcription factor，TF）。通过transmir（https://www.cuilab.cn）数据库对TF和miRNA之间的相互作用进行验证。最后通过cytoscape对TF，miRNA和hub基因之间的关系进行结果可视化。
1.6 Hub基因验证
从GEO数据库（https://www.ncbi.nim.nih.gov/geo）中下载基因微阵列数据集GSE206848，对Hub基因的表达水平进行验证。GSE206848基于GPL570[HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array检测平台，包含10例正常滑膜组织和7例OA滑膜组织。采用t检验分析两组间Hub基因表达水平的差异，并用受试者工作曲线分析Hub基因诊断OA的价值。以P＜0.05为差异有统计学意义。
1.7 免疫浸润分析
利用ImmuCellAI（http://bioinfo.life.hust.edu.cn）算法对GSE55235数据进行分析，根据P＜0.05过滤样本并计算每种免疫细胞的百分比。使用SangerBox生信分析绘制免疫细胞组成的柱状图及22种免疫细胞浸润可视化箱图，并分析正常组和OA组22种免疫细胞浸润的差异；使用仙桃学术在线绘图工具绘制22种免疫细胞相关性热图；使用微生信在线绘图工具绘制22种免疫细胞与OA组hub基因相关性热图。
2 结果
2.1 骨关节炎糖酵解相关差异基因的识别/筛选
GSE55235数据集分析共发现正常对照组和OA组之间共有444个差异基因，包括上调基因222个，下调基因222个（图1（a））。186个糖酵解相关基因从GSEA在线数据库获得，通过韦恩图对差异基因和糖酵解相关数据库取交集得到14个糖酵解相关差异基因（图1（b））。14个糖酵解相关基因热图和相关性热图分别见图1（c）和图1（d）。
图1骨关节糖酵解相关差异基因筛选
Fig.1Genetic variations associated with joint glycolysis screening
注：（a）GSE55235数据集的DEGs;（b）火山图糖酵解相关差异基因韦恩图; （c）糖酵解相关基因热图;（d）糖酵解相关基因相关性热图.
2.2 GO和KEGG富集分析
糖酵解相关差异基因筛选后富集到30个GO条目（P＜0.05），其中生物学过程10个，细胞组成10个，分子功能10个。应用R绘制CC，MF和BP气泡图，见图2（a）~2（c）。GO富集生物学过程主要包括细胞迁移、细胞活性、细胞定位、运动、细胞或亚细胞成分的运动、对缺氧的反应、对氧气水平的反应、磷脂酶a2活性的激活、磷脂酶a2活性的正调控; 细胞组成主要包括胞外区部分、细胞外基质、含有胶原蛋白的细胞外基质、基底膜、包膜囊泡膜、胶原V型三聚体、微脉冲、多囊体、内囊、纤维胶原三聚体、带状胶原纤维；分子功能主要包括肝素结合、糖胺聚糖绑定、硫化合物结合、14-3-3蛋白结合、电子转移活度、苹果酸脱氢酶（NAD+）活性、KDEL序列结合、苹果酸酶活性、草酰乙酸脱羧酶活性、谷胱甘肽二硫氧化还原酶活性等。同时，KEGG富集分析得到68条信号通路（P＜0.05），富集前10信号通路的弦图如图2（d）。与OA相关的信号通路主要有PI3K-Akt信号通路、EGFR酪氨酸激酶抑制剂耐药性、PI3K-Akt信号通路、自噬和HIF-1信号通路等。糖酵解相关基因GO和KEGG富集图形见图2。
图2糖酵解相关差异基因GO和KEGG分析（P＜0.05）
Fig.2Analysis of glycolytic related differential genes GO and KEGG（P＜0.05）
注：（a）GO（CC）富集分析结果（P＜0.05）；（b）GO（MF）富集分析结果（P＜0.05）；（c）GO（BP）富集分析结果（P＜0.05）；（d）KEGG富集分析结果（P＜0.05）.
2.3 PPI网络和核心基因筛选
采用STRING在线数据库构建了糖酵解相关差异基因的PPI，见图3（a），包括14个蛋白节点和9次相互作用，平均度值为1.29。利用Cytoscape软件中的cytohubba插件采用Degree算法筛选出连接度最高的8个基因，包括STC1，SDC3，IRS2，EGLN3，DDIT4，VEGFA，EGFR和ANG（图3（b））。
2.4 TF-miRNA-hub genes 调控网络构建
将8个Hub基因输入miRWalk数据库预测相应的miRNA，再通过mirTarBase，Targetscan，和miRDB数据库进行取交集筛选，获得30个miRNA，见图4（a）。通过TRRUST数据库筛选Hub基因相应的TF，共获得21个TF，见图4（b）。形成以EGFR，VEGFA，EGLN3，DDIT，SP1与STAT3为核心的互作网络，见图4（c）。
2.5 Hub基因验证
选择GSE206848数据集对Hub基因进行验证，结果显示hub基因在数据集中的表达趋势与GSE55235相一致，见图5（a）；且hub基因均具有较好的预测功能，几乎所有hub基因去曲线下面积均＞0.7，见图5（b）。
图3糖酵解相关差异基因PPI网络构建与模块分析
Fig.3PPI network construction and module analysis of differential genes related to glycolysis
注：（a）STRING在线数据库构建GDEGs的PPI；（b）Degree算法筛选的连接度最高的基因.
图4TF-miRNA-hub genes 调控网络构建
Fig.4TF-micrornas-the hub genes regulating network building
注：（a）Hub-miRNA网络；（b）Hub-TF网络；（c）Hub-TF-miRNA三层网络.
图5糖酵解相关差异基因Hub基因验证
Fig.5Validation of Hub gene related to glycolysis
注：（a）GSE206848数据集Hub基因箱图；（b）GSE206848数据集Hub基因ROC图.
2.6 免疫浸润分析
22 种免疫细胞在样本中的百分比展示在图6（a）。与正常滑膜相比，浆细胞和静息型肥大细胞在OA组织中浸润深度较高，差异有统计学意义（P＜0.05）；记忆CD4T细胞、调节性T细胞、激活型肥大细胞和嗜酸粒细胞在OA组中浸润深度显著低于对照组，差异有统计学意义（P＜0.05），见图6（b）。EGFR，VEGFA，DDIT4，EGLN3，IRS2和STC1与静息型T细胞以及激活型肥大细胞呈正相关，与静息型肥大细胞呈负相关；SDC3与静息型T细胞、巨噬细胞M1以及激活型肥大细胞呈负相关，与浆细胞、巨噬细胞M2及静息型肥大细胞呈正相关，见图6（c），6（d）。
图6免疫浸润分析图
Fig.6Immunoinfiltration analysis diagram
注：（a）GSE55235数据集22种免疫细胞在样本中的百分比柱状图（P＜0.05）；（b）GSE55235数据集22种免疫细胞表达趋势可视化箱图（P＜0.05）；（c）GSE55235数据集22种免疫细胞相关性热图；（d）GSE55235数据集22种免疫细胞与hub基因相关性热图.
3 讨论
以前学者对于OA的研究大多数集中于关节软骨，而膝关节周围软组织的功能却被忽略了[8]。近期越来越多的研究发现，OA的早期至晚期均有可能发生滑膜炎[9]，然而目前对OA滑膜炎了解远远不够，因此，研究滑膜炎的机制也至关重要。此外，尽管在很多疾病中糖酵解已经得到了广泛的研究，但在滑膜炎中，其功能和作用仍然不太清楚。
本研究将健康者与OA 病例的滑膜样本相比较后的差异基因与糖酵解相关基因相交获得14个糖酵解相关差异基因，之后进行富集分析。这些糖酵解相关差异表达基因（Glycolysis related differentially expressed genes，GDEGs）显著富集了缺氧反应与酶活性，表明这些GDEGs与氧化有关，与之前的研究结果一致[10]。另外，通过KEGG通路富集，我们发现这些基因主要富集于HIF-1信号通路和PI3K-AKT中。HIF-1α是调节缺氧发育和细胞反应的主要转录因子。已有研究证明HIF-1α的丢失会加重MMP13的水平并降解软骨组织[11]。
PI3K/AKT 信号通路在许多细胞存活通路中起着重要的调节作用，主要作为细胞凋亡的抑制剂。研究发现PI3K/Akt 信号通路参与许多生物过程，包括细胞周期、细胞凋亡、血管生成和葡萄糖代谢。此外，PI3K/Akt还可以减少糖原的合成并增加糖酵解[12]。而功能失调的 AKT 网络与多种病理状况有关，包括过度生长综合征、癌变、心脏病、糖尿病、自身免疫性疾病和炎症以及神经系统疾病。现有研究表明PI3Kδ和PI3Kγ在类风湿关节炎滑膜和培养的滑膜细胞中高表达，与滑膜炎症的调节密切相关[13]。这或许可以帮助我们增加对OA机制的了解。
然后，使用Cytoscape中的cytoHubba程序检验了PPI结果，得到了8个关键基因，即STC1，SDC3，IRS2，EGLN3，DDIT4，VEGFA，EGFR和ANG。其中EGFR已经被证实对软骨细胞有双重调节作用。其主要在出生后和成人关节软骨中起合成代谢作用。在OA损伤时，需要EGFR信号传导来限制软骨退化，并可能用于再生[14]。VEGFA作为枢纽基因富含生长因子受体结合和生长因子活性，它通过刺激血管内皮细胞的有丝分裂和诱导内皮细胞的迁移来促进血管生成。它可以增强和维持毛细血管的通透性，不仅有利于萌芽形式的新生血管的形成，而且有利于营养物质的运输和炎症的扩散[[15]。适当的VEGFA水平可以改善滑膜细胞对关节腔缺氧环境引起的细胞损伤的适应。然而，VEGFA的过度表达会导致滑膜的异常增殖和滑膜炎的加重[16]。这些关键基因经过验证GSE206848表明这些基因对OA也具有较可靠的诊断价值。除此之外，这些关键基因经过KEGG分析主要富集于HIF-1信号通路和PI3K-AKT中。有研究表明，HIF-1α可以激活相关糖酵解酶的表达，然后通过PI3K/Akt途径调节糖酵解[17]。但是HIF-1α调控的糖酵解靶基因与PI3K/Akt调控的关联仍有待探索。
富集通路分析还显示miRNA 也与 OA 有关，并且几种 miRNA 与疾病发病机制有关[18]。因此本研究建立了miRNA靶基因网络，在此基础上，发现hsa-miR-7-5p和hsa-miR-200a-3p与验证后最具有诊断预测性的EGFR相互作用。hsa-miR-7-5p在类风湿性关节炎（RA）患者的滑膜中表达水平升高[19]，因为RA是一种慢性全身性炎症性自身免疫性疾病，伴有骨和关节破坏，所以大部分研究认为其与骨关节炎具有密切联系，但hsa-miR-7-5p在RA中的潜在机制还需要进一步阐明。hsa-miR-200a-3p在炎症环境中表达显著，但目前没有确切的研究表明它与OA的机制有关联。根据本文的研究结果，这些miRNA可能参与滑膜的糖酵解，可它们在OA滑膜增生中的作用尚不明确，需要更多的研究来进一步阐明。
另外据研究表明，糖酵解与类风湿性关节炎等免疫疾病也有关联[20]，但是现在糖酵解和免疫之间的确切机制还不清楚，因此，使用CIBERSORTx来研究OA滑膜中免疫细胞的浸润。结果显示，OA滑膜中的一些免疫细胞与正常滑膜有显著差异。与正常滑膜相比，OA滑膜中的T调节细胞（Tregs）、巨噬细胞M0和静息型肥大细胞的比例相对较高，而活化CD4记忆T细胞和活化肥大细胞的比例较低[21]。巨噬细胞在损伤后伤口愈合反应中作为先天免疫和组织重塑的关键细胞介质已经确定[22]，在OA中巨噬细胞似乎也起着类似的作用。很多研究证明滑膜巨噬细胞与软骨降解、软骨变性以及骨赘生成有着密切的联系[23]。在Wood等[24]。的研究中发现OA患者滑膜中的免疫细胞绝大多数由巨噬细胞主导，而T细胞在炎性关节炎患者的滑膜中占主导地位，这与本研究的免疫浸润结果一致，也证明了炎症表型或许是OA治疗的一种新思路。另外，几项临床研究报告了OA患者滑膜组织中参与肥大细胞分化和活性的基因表达增加，保护MC缺陷小鼠免受OA的炎症和软骨破坏[25]，表明了肥大细胞可能也是OA的潜在治疗靶点[26]。
4 结论
综上所述，本研究采用生物信息学方法对数据集进行分析，在正常滑膜组织与OA患者滑膜组织之间筛出DEGs，进一步筛选出8个关键基因并进行验证，最终筛选出EGFR，VEGFA，EGLN3和DDIT4四个糖酵解表型差异基因为OA期间滑膜炎潜在的诊断生物标志物和治疗靶点，为临床诊断和治疗提供参考。但这项研究仍然存在一些局限性。首先，研究使用的数据集可能会存在偏差，这将会导致最后的结果也随之发生偏倚。其次，我们还需要更多的临床研究或实验来验证这些因素是否有诊断和预后意义。最后，全面了解OA的关键基因靶点对于完善OA的发病机制和促进制定相对完善的治疗策略是至关重要的。
