骨关节炎相关基因差异的表达分析

陈 璞

（东莞市第八人民医院骨关节科，广东 东莞 523000）

骨性关节炎（osteoarthritis，OA）是世界范围内的慢性疾病，以前骨关节炎只是简单的被认为由于磨损而受到的损害，而现在一致认为骨关节炎是由局部和系统因素的复杂作用造成的［1］。骨性关节炎发生的年龄一般为40 岁以上，主要表现在手指、肩膀、臀部、膝盖或脚踝疼痛等。50 岁以上的患者有关节疼痛、晨起僵硬和功能障碍等表现［2］。这些骨性关节炎的发生严重影响老年人的身心健康，给患者和政府医疗带来了巨大的经济负担。

骨性关节炎的发病与多种因素有关，关节滑膜对关节软骨的生长有一定的促进作用，有研究表明OA 发生的重要机制是由滑膜病变引起的［3］。在OA发生与发展当中，滑膜病变通过多种途径引起骨类疾病的发生与发展［4］。近年来OA 的遗传学研究成果丰硕，表观遗传学、基因组学和蛋白质组学研究表明，OA 相关的基因参与了滑膜发育相关通路，这些基因的突变导致有关通路的紊乱最终可能会造成OA 的发生［5］。目前，虽然有大量研究正在探索OA 的发生机制，但是其真正的发病原因及机制尚不明确，其发病机制的研究对OA 的早期诊断及有效的治疗具有重要的意义。目前主要采用影像学检查来进行早期的诊断，由于早期诊断的手段比较单一，因此，骨性关节炎的分子机制研究显得尤为重要。

随着高通量测序技术的发展，生物学领域的公共数据库共享了大部分与疾病相关的基因表达谱数据，这些海量数据很容易获得，包括本文研究的OA 数 据 集。 本 文 利 用R 语 言“limma”包 对GSE55235［6］和GSE55457［6］芯片预 处理、差异 分析等，筛选出OA 患者与正常患者的上调和下调的差异基因，接着对差异表达基因做GO和KEGG富集分析，然后通过R 语言“STRINGdb”包构建了基于差异表达基因的PPI网络图，以此筛选与OA 发生与发展过程中的关键基因。本研究的主要目的是筛选与OA 的发病机制相关的关键基因，并揭示OA 的潜在的发病机理。

1 材料与方法

1.1 GEO 芯片数据的获取从美国国立卫生研究院的NCBI 下载基因表达数据库（Gene Expression Omnibus，GEO），GEO（https：//www.ncbi.nlm.nih.gov/geo/）的基因表达谱是现在最全面的公开数据集［7］。从NCBI 下载的基因表达谱芯片数据GSE55235 与GSE55457，这两套芯片均为Affymetrix human genome U133A 芯片，平台都为Affymetrix GPL96，2 个数据集都包含来自正常膝关节的10 例滑膜组织样本以及关节置换下来的10 例膝关节骨性关节炎滑膜组织样本。

1.2 差异表达基因的筛选芯片数据经过预处理后进行差异表达分析。利用R 语言“affy”［8］包进行预处理和“limma”［9］包进行差异表达分析，得出差异表达基因。其中调整后的P 值（adjusted p）＜0.05，|logFC|＞1 作为差异基因的筛选条件分别选择两个数据集中差异表达的基因，并且将两个数据集中表达上调或下调的基因分别取交集筛选出共同表达上调或下调的基因作为进一步分析的对象，得到了差异基因的表达热图。我们采用多重检验校正（FDR）来校正P值。

1.3 差异基因的GO 和KEGG 富集分析采用R语言“clusterProfiler”［10］包对与骨性关节炎相关的差异基因进行了GO 和KEGG 富集分析。其中GO 富集分析分别从生物进程（Biological processes，BP）、细 胞 组 成（Cellular component，CC）和 分 子 功 能（Molecular function，MF）等3 个 方 面 进 行 富 集。KEGG 富集分析差异基因相关通路，计算这些差异基因与通路的相关显著性，筛选到与OA 发生与发展密切相关的信号通路。通路以基因数目＞2 和富集的P.adjust＜0.05作为富集分析的筛选标准。

1.4 PPI 网络构建STRING 数据（http：//www.string-db.org/）是蛋白质相互作用网络的数据库［11］。其包含蛋白质之间的直接相互作用和间接相互作用。本研究通过R 语言“STRINGdb”［12］包构建了骨性关节炎差异基因的PPI网络。

2 结 果

2.1 差异基因分析开始使用affy 包进行原始芯片归一化处理，然后通过R语言“limma”包的差异分析，GSE55235 数据集共筛选出差异基因总数为756，其中包括467 个上调基因和289 个下调基因。GSE55457 数据集共筛选出240 个差异基因，其中上调基因56 个，下调基因184 个。将两个数据集中得到的上调和下调的基因分别取交集，共得到共同上调基因20 个，共同下调基因36 个。基于共有的56个差异基因，表1 中显示了前30 个差异表达基因的表达情况，这些差异基因是以adjusted P＜0.05 和|logFC|＞1为标准筛选的。利用R语言中的“pheatmap”包绘制GSE55235 和GSE55457 基因热图（图1～图2）。

表1 前30个差异表达基因的显著性情况

图1 GSE55235差异表达基因热图

图2 GSE55457差异表达基因热图

2.2 GO 富集分析为了探索差异基因的生物学功能，我们通过R 语言clusterProfiler 包对骨性关节炎差异基因进行了GO 富集分析，由BP、CC 和MF三部分构成（图3）。生物进程（BP）中差异基因主要富集在蛋白质磷酸化的负调控（negative regulation of protein phosphorylation）、磷酸化的负调控（negative regulation of phosphorylation）、上皮细胞增殖（epithelial cell proliferation）、细胞对外部刺激的反应（cellular response to external stimulus）、转移酶活性的负调控（negative regulation of transferase activity）中。在细胞组成（CC）中，主要富集在蛋白质细胞外基质（pro‑teinaceous extracellular matrix）、细胞外基质（extra‑celluar matrix）中，关节软骨是透明软骨，由细胞外基质包裹软骨细胞组成，细胞外基质的主要成分是胶原和蛋白多糖，正常软骨细胞外基质合成和降解处于动态平衡，当基因表达异常的时候，平衡状态被破坏就会导致OA 的发生。在分子功能（MF）中，主要富集在糖胺聚糖结合（glycosaminoglycan binding）、生长因子结合（growth factor binding）、MAP 激酶磷酸酶活性（MAP kinase phosphatase activity）、糖皮质激素受体结合（glucocorticoid receptor binding）中。

2.3 KEGG 富集分析为了筛选到与OA 相关通路，通过clusterProfiler 包对骨性关节炎相关的差异基因进行了KEGG 通路富集分析（图3）。通过KEGG 富集分析，发现差异基因主要集中在MAPK signaling pathway、Epstein-Bar virus infection、Human T-cell leukemia virus 1 infection、PI3K-Akt signaling pathway、 Kaposi sarcoma-associated herpesvirus infection、ErbB signaling pathway 和HIF-1 signaling pathway 等通路中。WANG 等研究表明乳铁蛋白通过MAPK信号通路促进骨关节炎老鼠的软骨细胞增殖［13］。ZHANG 等证实microRNA‑130a 通过PTEN/PI3K/Akt信号通路引起软骨细胞骨关节炎［14］

图3 OA差异基因GO功能分析和KEGG通路分析以及差异基因构建的PPI网络

2.4 差异基因的PPI 网络分析基于STRING 的蛋白质网络数据库，我们通过R 语言“STRINGdb”包构建了PPI 网络，其中包括39 点，95 个边（图3D）。在56 个差异基因中，进行网络拓扑学分析，按照度的大小来表示该基因在网络中的重要性的标准，度值排在前10 的基因分别为VEGFA、JUN、MYC、ATF3、DUSP1、NR4A1、CDKN1A、NR4A2、GADD45B和KLF6（表2）。

表2 PPI网络中前10基因的度情况

3 讨 论

OA 是骨科常见的关节慢性疾病，在美国的发病数目达到千万以上［15］。其特点主要是滑膜无菌性炎症，其临床症状为关节肿胀、疼痛等。甚至有可能导致残疾，使得中老年人的生活质量受到极大的影响，并且给患者带来巨大的经济压力［16］。所以对患者的早期诊断与治疗非常重要。

本研究通过对OA 表达谱芯片GSE55235 和GSE55457进行生物信息学的差异表达分析，总共得到56 个差异基因，其中包括20 个表达上调基因和36 个表达下调基因。然后对得出的差异表达基因进行GO 和KEGG 富集分析，得到与OA 发病相关的通路。最后，通过R 语言“STRINGdb”包来构建了差异表达基因的PPI 网络图，其中包括39 个基因和95条相互关系。在56 个差异基因中，筛选出VEGFA、JUN、MYC、ATF3、DUSP1、NR4A1、CDKN1A、NR4A2、GADD45B和KLF6基因作为该网络的关键基因。有报道表明PI3K-Akt 与MAPK 信号通路都属于破骨分化信号通路［17］，c-JUN（JUN）氨基末端激酶（c-JUN N-terminal kinase，JNK）作为MAPK 信号通路的关键分支，成为炎症药物治疗的常见靶点［18］。而且在白细胞介素1 的诱导下JNX/c-JUN 信号通路可以抑制Col-2基因的表达，然后关节软骨细胞外基质正常组分的合成受到影响，从而降低关节软骨的强度及韧性，最终导致骨关节炎的发生［19］。FISCH 等选择了丰富的转录因子来调节OA 中的差异表达基因，并通过创建软骨特异性相互作用的局部网络来对这些转录因子进行了优先排序。该分析揭示了8个转录因子，其中包括MYC 作为OA 中大量失调基因的靶标之一，并在OA中被抑制表达［20］。Atf3通过调节软骨细胞中的炎性细胞因子的表达来参与OA 的发生与发展，而软骨细胞中的炎性细胞因子/NF-kB/Atf3 的前馈环可能是OA 治疗的新型治疗靶点［21］。有研究表明与正常细胞相比，CDKN1A 基因（细胞周期蛋白依赖性激酶抑制剂1A，也称为p21 或野生型p53 激活片段1）在OA 软骨细胞中被下调［22］。由于其能够阻止许多细胞类型增殖，因此提出了CDKN1A 的缺失导致OA 软骨中细胞增殖的重新启动［23］。目前为止并没有发现KLF6 与OA 相关联，但是实验发现KLF6 基因在骨关节炎与正常组织对比中表达具有显著的差异，因此KLF6 可能在骨关节炎的发生发展中起到一定的作用，其具体的作用机制有待进一步研究。血管内皮生长因子（VEGFA）水平升高与骨关节炎（OA）的发展有关，实际上VEGF可能参与了OA 特有的病理进程，包括软骨变性、骨赘形成、软骨下骨囊肿和硬化、滑膜炎和疼痛。此外，广泛的研究表明抑制VEGF信号传导可降低OA的发展［24］。针对VEGFA、JUN、MYC、ATF3、DUSP1、NR4A1、CDKN1A、NR4A2、GADD45B和KLF6基因的检测将有助于OA 的早期诊断与治疗。这也预示着这些基因可能与OA的发病机制有关。

通过生物信息分析方法筛选到了OA 的差异表达基因，并对其进行了富集分析，然后通过PPI网络分析，得出与OA 发生发展相关的相关基因。这有助于骨性关节炎的早期诊断和治疗，为寻找与OA发生发展及治疗相关靶点的研究提供了新的思路。