骨关节炎是一种以软骨变性、滑膜炎以及骨赘形成等为特征的慢性炎症和退行性关节疾病，常导致老年人慢性残疾。目前骨关节炎的致病机制尚不明确。以往骨关节炎的研究热点大多集中于软骨组织，随着研究不断深入，滑膜病变在骨关节炎发生发展中的作用引起人们的重视。研究发现超过90%的骨关节炎患者被证实存在滑膜病变，且骨关节炎滑膜的病变程度与关节严重的疼痛和功能障碍有关[1]。此外，滑膜炎的发生还可能促进软骨退变[2]。寻找滑膜病变的产生机制有助于找到精确的治疗靶点，获得良好预后，从而从根本上解决骨关节炎问题，这有赖于生物学标志物的鉴定。因此寻找疾病差异基因、相关信号通路是骨关节炎发病机制及治疗靶点研究中亟待解决的问题。

随着基因芯片和RNA测序技术的迅速发展，生物信息学分析成为一个重要的研究方向，为识别可靠的功能差异表达基因(differentlyexpressedgenes,DEGs)、微小RNA(miRNAs)、环状RNA(cirRNAs)和长链非编码RNA(IncRNAs)提供了新的线索和核心数据，在筛选各种疾病候选生物标志物方面也展示出了巨大的优势。然而，由于样本源于不同的测序平台，表达的mRNA结果与基因谱对应不一致，很大一部分骨关节炎的生物信息学分析仅限于单一芯片数据，结果较为局限，可靠性差，而多组基因芯片数据可有效地解决这一问题。此外，随着大量研究的进行，大量上传至公共数据库的基因信息未被有效利用。因此，此次研究利用多组基因芯片数据再次进行筛选与分析，拟鉴定出更具有生物学意义的标志物，为骨关节炎的研究提供新思路。

1、资料和方法

1.1 资料来源

研究从公共数据库GEO(https://www.ncbi.nlm.nih.gov/geo/)下载骨关节炎滑膜相关芯片数据集GSE82107、GSE12021、GSE55457、GSE55235，筛选条件为：(1)骨关节炎；(2)人类；(3)滑膜；(4)缺乏药物干预。其中纳入GSE82107数据集10例骨关节炎滑膜组织，7例正常滑膜组织；GSE12021数据集10例骨关节炎滑膜组织，9例正常滑膜组织；GSE55235数据集10例骨关节炎滑膜组织，10例正常滑膜组织；GSE55457数据集10例骨关节炎滑膜组织，10例正常滑膜组织。共计纳入骨关节炎滑膜组织40例，正常滑膜组织36例，其芯片平台包括GPL570[HG-U133＿Plus＿2]AffymetrixHumanGenomeU133Plus2.0Array和GPL96[HG-U133A]AffymetrixHumanGenomeU133AArray两种类型。

1.2 分析方法

1.2.1 差异表达基因筛选

用R语言软件读取下载的原始数据，采用Oligo等软件包进行数据处理，sva包去除4组数据间的批次效应，获得经过批次校正后的新数据矩阵，并标注矩阵中数据的骨关节炎组和健康对照组。采用Limma包进行差异表达分析，设定P<0.05,adjustPValue<0.05,log2基因表达差异倍数(Fold-change,FC)绝对值≥0.58为筛选条件，筛选获得差异表达基因。

1.2.2 GO和KEGG功能富集分析

利用在线数据库对差异表达基因进行功能和通路富集分析，进一步探讨疾病的发病机制。将筛选获得的差异表达基因载入至R软件中，使用org.Hs.eg.db、ggplot2、clusterProfiler软件包，enrichGO、enrichKEGG算法进行基因本体论(Geneontology,GO)和KEGG(Kyo-toEncyclopediaofGenesandGenomes)信号通路分析。

1.2.3 PPI网络构建及关键基因筛选

利用在线分析网站寻找核心基因。将所有的差异基因差异表达基因上传至STRING在线分析网站，导出最低可信度“combinedscore>0.7”的分析结果。之后进一步通过Cytoscape得到差异基因的蛋白-蛋白交互作用网络模型(protein-proteininteractionnetwork,PPInetwork)，同时采用其中的Hubba插件，采用Degree、Betweenness、Closeness三种不同算法，分别计算出各自得分排名靠前的20位基因，并用在线分析平台DrawVennDiagram(http://bioinformatics.psb.ugent.be/webtools/Venn/)将三者取交集，得出关键基因。

2、结果

2.1 差异表达基因

在进行差异基因的分析之前，需要先对单个芯片数据预处理。如存在所有样本中都不表达的探针、缺失值存在或基因与多个探针存在对应关系等特殊情况，需进行标准化、探针过滤、缺失值填充以及探针合并等。随后考虑到4组数据芯片间批间差异的存在，去除批次效应。主成分分析PCA(PrincipalComponentAnalysis)结果可以看出，批间差处理之后的芯片数据分布效果较好(图1)。

进一步进行差异分析，结果显示对比骨关节炎组与健康对照组筛选得到447个差异基因，其中骨关节炎组中高表达基因201个，低表达基因246个。对差异基因进行可视化分析绘制火山图和热图。火山图中每个点代表一个基因，其中基于P值<0.05和|logFC|≥0.58的标准，蓝色表示低表达的基因，红色表示高表达的基因(图2)。热图直观展示了骨关节炎与健康对照组之间的基因表达差异(图3)。骨关节炎患者滑膜差异表达基因见表1。

2.2 GO功能富集分析

GO分析主要包括生物学过程、细胞组成以及分子功能3个部分。其中上调的差异表达基因主要富集在炎症反应、细胞外基质、胶原蛋白、胶原蛋白三聚体、配受体活性等方面；下调的差异表达基因主要富集在血管发育、脂质代谢过程的调控等生物学过程(图4)。

2.3 KEGG信号通路分析

KEGG结果显示上调的差异表达基因主要集中在吞噬体、类风湿性关节炎、产IgA的肠免疫网络、金葡菌感染、同种异体移植排斥反应、花生四烯酸代谢等信号通路；下调的差异表达基因主要集中在MAPK信号通路、肿瘤坏死因子、白细胞介素17、低氧诱导因子1信号通路、癌症转录失调等信号通路(图5)。

2.4 差异表达基因的PPI分析

基于Cytoscape软件，筛选出degree、betweenness、closeness得分高的前20位基因，将三者取交集筛选出基质金属蛋白酶9(MMP9)、白细胞介素6(IL-6)、血管内皮生长因子A(VEGFA)、JUN、前列腺素内过氧化物合酶2(PTGS2)、CXCL8、MYC、表皮生长因子受体(EGFR)8个关键基因(图6)。KEGG功能富集分析显示，这些关键基因主要在白细胞介素17信号通路、肿瘤坏死因子信号通路、巨细胞病毒感染和癌症相关通路富集(图7)，同时表2展示了这些通路分别对应的关键基因。

3、讨论

此次研究从GEO数据库提取滑膜原始数据，使用生物信息学分析挖掘出骨关节炎与健康对照之间的差异表达基因，其中上调基因有201个，下调基因有246个。对这些差异基因进行功能富集和通路分析以及蛋白相互作用分析之后，筛选出基质金属蛋白酶9、白细胞介素6、血管内皮生长因子A、JUN、PTGS2、CXCL8、MYC、表皮生长因子受体8个关键基因。

GO分析显示下调的差异表达基因主要富集在血管发育方面。其中血管内皮生长因子A是骨发育过程中的重要调节因子[3]，骨关节炎晚期的关节软骨和滑膜中血管内皮生长因子表达明显增加。它参与了骨关节炎的软骨变性、骨赘形成、软骨下骨囊肿和硬化、滑膜炎等特异性病变过程，在白细胞介素1β的协同下,血管内皮生长因子能够显著降低聚集蛋白聚糖和Ⅱ型胶原蛋白的基因表达和蛋白质水平[4]，而抑制血管内皮生长因子A信号传导可延缓骨关节炎进展。有意思的是，最新研究提示血管内皮生长因子A信号通路的功能存在两面性，它既可以促进血管新生也可以诱导血管退化[5]。此次结果显示血管内皮生长因子A基因在骨关节炎滑膜中是低表达的，与以往的研究不一致，作者猜想是否由于机体内血管内皮生长因子蛋白的过表达反向抑制血管内皮生长因子A基因的表达，从而导致基因下调，同时也不能排除结果的不一致可能与基因库的选择不同有关。总之，血管内皮生长因子A在骨关节炎的发生发展中起着重要的作用，而其发挥作用的具体机制需要进一步进行探究。

表1|骨关节炎患者滑膜差异表达基因

表2|关键基因KEGG功能富集分析

基质金属蛋白酶9是一种明胶酶，它参与正常生理过程中细胞外基质的代谢过程，如胚胎发育、生殖和组织重塑，也参与一些疾病如关节炎和癌症的发生过程。循证医学证据Meta分析结果提示骨关节炎患者基质金属蛋白酶1、基质金属蛋白酶2、基质金属蛋白酶9蛋白水平高于健康对照组，并且基质金属蛋白酶9蛋白水平的升高不受种族的影响[6]。许多研究均已证实了基质金属蛋白酶9参与骨关节炎的发病过程[7]，基质金属蛋白酶9对Ⅳ型胶原蛋白具有蛋白水解活性，通过降解骨关节炎患者关节内非胶原基质成分，使关节内结缔组织遭到破坏。与以往大多数研究结果一致，此次研究发现骨关节炎患者的基质金属蛋白酶9基因显著上调，由此作者认为基质金属蛋白酶9很可能是骨关节炎重要的生物学标志物及潜在治疗靶点。此外还有学者发现基质金属蛋白酶9与血管内皮生长因子也存在密切的关系[8]，因此明确血管内皮生长因子A与基质金属蛋白酶9之间的关系，对探索骨关节炎的治疗药物具有重要作用。

图1|批间差处理前后PCA图

图2|基因表达的火山图

图4|差异基因GO功能富集分析气泡图

图3|前50名上调和下调差异表达基因的热图

图5|差异基因KEGG功能富集分析气泡图

图6|cytoscape软件3种不同算法下的Top20

图7|关键基因KEGG功能富集分析气泡图

近年来，破骨细胞在慢性关节炎中的作用逐渐引起重视，破骨细胞可能是导致关节破坏和骨吸收的关键因素[9]。肿瘤坏死因子信号通路与MAPK信号通路均参与破骨细胞分化，促炎性细胞因子如肿瘤坏死因子α、白细胞介素1和白细胞介素17可以促进破骨细胞生成[10]。c-JUN氨基端激酶(c-JunN-terminalkinase,JNK)是MAPK信号转导通路的重要组成之一，其在骨关节炎关节软骨细胞中所起的作用已被证实[11]。c-JUN与碱性亮氨酸拉链转录因子ATF样蛋白(basicleucinezippertranscriptionfactor,ATF-like,BATF)结合形成的异二聚体复合物可以调节软骨细胞内合成和分解代谢基因的表达，参与骨关节炎软骨的破坏过程中起重要作用[12]。此外，MAPK通路中JUK的激活能够促进白细胞介素6和肿瘤坏死因子α的表达[13]，从而间接通过肿瘤坏死因子信号通路诱导破骨细胞分化。然而，此次研究中JUN基因显示下调与以往研究尚不统一，需进一步验证。

以往有研究显示骨关节炎患者滑膜成纤维细胞和滑膜液中白细胞介素6均过表达[14,15]。此次研究发现白细胞介素6基因在骨关节炎与健康对照组滑膜对比中表达减低。这可能与携带白细胞介素6的C等位基因个体对骨关节炎的易感性较低且损伤较轻，导致骨关节炎患者中白细胞介素6的表达降低有关[16]。也有人对血清及滑液中白细胞介素6表达水平与疾病的相关性进行研究。如RADOJČIĆ等[17]的研究结果显示，关节滑液中的白细胞介素6表达水平与患者的疼痛程度显著相关，而血清中白细胞介素6的表达水平与疼痛的严重度之间无明显的相关性；SHIMURA等[18]的研究也证实了这一点。然而对于骨关节炎患者滑膜中白细胞介素6表达量的高低以及是否与疾病进展相关联尚未定论，仍需进一步探讨。

此次研究KEGG分析结果显示，前列腺素内过氧化物合酶2参与的花生四烯酸代谢途径富集于差异表达基因的KEGG是上调的，提示骨关节炎患者滑膜组织中前列腺素2表达也升高。以往研究认为前列腺素E2是引起关节炎症性疼痛的主要因素[19]。作为一种促炎因子，环氧合酶2在白细胞介素1β和肿瘤坏死因子α的刺激下表达增加，诱导了前列腺素内过氧化物合酶2的产生，进而促进骨的吸收[20]。目前临床上使用的非类固醇抗炎药的作用机制主要是通过抑制环氧合酶来抑制致炎前列腺素的生物合成，减轻炎症反应，有良好的应用价值。那么抑制前列腺素内过氧化物合酶2是否也能够减轻患者病痛需要进一步探究。

此次研究也证实骨关节炎患者滑膜液中趋化因子CXCL8表达明显增加。趋化因子CXCL8，即白细胞介素8，是巨噬细胞和上皮细胞等分泌的细胞因子，其主要生物学活性是吸引和激活中性粒细胞从而发挥促炎反应。中性粒细胞与CXCL8接触后发生形态变化，定向游走到炎症反应部位并释放一系列活性产物；这些产物可导致机体局部炎症反应的发生，达到杀菌和破坏病菌感染细胞的目的。此外CXCL8也可能通过调控JAK-STAT、NF-kB和JNK-MAPK信号通路，促进其他促炎细胞因子的释放，诱导细胞凋亡，抑制软骨细胞的增殖，从而诱导骨关节炎的发生[21]。

表皮生长因子受体是表皮生长因子受体家族成员之一。表皮生长因子受体位于细胞膜表面，与配体表皮生长因子和转化生长因子α结合后激活可产生自磷酸化由单体转化为二聚体，其自磷酸化后可以引导下游信号通路(包括MAPK,Akt和JNK通路)的磷酸化，导致细胞增殖。表皮生长因子受体也可以通过调节血管生成素1及血管内皮生长因子等因子表达水平而影响血管生成，虽然大多数数据都支持表皮生长因子受体信号途径在软骨内成骨和关节软骨生长中起重要作用[22,23]。然而，也有研究认为表皮生长因子受体激活的可能仅仅是骨关节炎的某一细胞亚群，约占搜集样本的27%[24]，此次结果也显示表皮生长因子受体基因在骨关节炎滑膜中低表达。此外表皮生长因子受体对骨关节炎啮齿类动物模型的影响可能会有矛盾的结果产生[22,25]。因而针对此治疗靶点的研究仍需进一步探索。而对于MYC的探讨大多研究集中于癌症相关疾病，与骨关节炎的关系尚不清晰。

综上所述，此次研究揭示包括基质金属蛋白酶9、白细胞介素6、血管内皮生长因子A、JUN、前列腺素内过氧化物合酶2、CXCL8、MYC、表皮生长因子受体在内的hub基因在骨关节炎中发挥重要作用，而这些生物学标志有可能成为骨关节炎的药物靶点和诊断标志物，其相关通路也具有重要的研究价值。而且这些基因所表达的蛋白分子相互联系，考虑是否可以通过对这些靶点的共同检测使诊断更加准确，这将有待于进一步研究。研究为骨关节炎的分子生物学研究提供了新的认识，也为发现诊断骨关节炎的新的生物标志物铺平了道路。

宋珊,胡方媛,乔军,王佳,张升校,李小峰.基于生物信息学途径认识骨关节炎滑膜的生物学标志物[J].中国组织工程研究,2021,25(05):785-790.