基于生物信息学方法鉴别骨关节炎自噬的关键基因和途径

袁长深 容伟明 李 哲 官岩兵 梅其杰 段 戡

1.广西中医药大学第一附属医院骨三科，广西南宁 530023；2.广西中医药大学研究生院，广西南宁 530000

骨关节炎（osteoarthritis，OA）是一种伴随结构改变的关节疼痛综合征，是最常见的关节疾病。目前其发病机制不明确，缺乏有效手段逆转该病的进程[1]，其导致的关节疼痛、残疾问题较为突出，加之该病容易丧失劳动力，故将其防治作为区域卫生政策的优先事项对待，甚至建议通过修订有关OA 临床和公共政策加以重视[2]。

软骨细胞死亡/存活在OA 发病中起着关键作用。在OA 的软骨细胞中，自噬相关标志物的表达显著降低，表明自噬可能是OA 重要的发病机制[3]，甚至决定OA 的命运[4]。自噬虽不是细胞死亡的方式，但可通过清除受损和功能失调的细胞器和大分子，维持OA软骨内环境的稳定。通过增强体内自噬平衡机制，可作为延缓关节衰老和降低OA 风险的一种新途径[5]。

目前关于OA 自噬的研究仍处于早期阶段，利用生物信息学方法挖掘自噬参与OA 发生、发展的重要基因及通路，可从一个新角度更好地阐明OA 发病机制，也为开发新药提供可靠的理论依据。

1 材料和方法

1.1 数据来源

通过美国国立生物技术信息中心建立的GEO 数据库中的GEO DataSets 检索窗口，以“Osteoarthritis”为检索词，同时限制条目类型为DataSets、研究的类型为Expression profiling by array、种属类型为homo sapiens 和公开日期为2010 年1 月1 日至2020 年1 月1 日，基于GPL96 平台，最终筛选得到基因微阵列数据集（GSE55457）；该数据集包含10 例正常和10 例OA 关节的滑膜组织，将其分别作为对照组和实验组。

1.2 差异表达基因筛选

使用R 软件limma 包对20 例样品（对照组和实验组）的表达矩阵进行校正和分位数归一化，再利用高级的t 检验分析并根据假设检验对P 值进行校正，筛选出校正后P ＜0.05 和|logFC|>1 的基因作为差异表达基因（differentially expressed genes，DEGs），并构建火山图。

1.3 OA 自噬基因的筛选

通过GeneCards 网站使用“自噬”检索词获取自噬相关基因。由于相关性评分被用来表示基因与自噬活性之间的强度，本研究将相关性评分≥1 分的基因作为自噬相关基因。将差异分析获得的DEGs 与自噬相关基因进行交集，利用R 软件Venn Diagram 包绘制交集图，交集基因即为OA 自噬基因。

1.4 GO 和KEGG 富集分析

对筛选出的基因进行GO 和KEGG 富集分析，可清楚其在疾病中的主要功能和途径，有助于深入研究疾病的机制。通过R 软件中Bioconductor 软件包对交集基因进行GO 和KEGG 富集分析，筛选出校正后P ＜0.05 的结果，再利用ggplot2 包对显著富集的结果分别进行可视化。

1.5 蛋白质-蛋白质相互作用（proteinprotein interaction，PPI）网络构建与模块分析

通过STRING 数据库在线构建交集基因的PPI网络，并设定筛选阈值为0.4，获得交集基因编码的PPI网络图，再将PPI 网络数据导入Cytoscape 软件，并使用MCODE 插件进行网络聚类，以筛选出主要模块，最后对其基因再次做KEGG 富集分析，比较两种基因所得信号通路的异同。

1.6 Hub 基因的选择

使用CytoHubba 插件对PPI 网络进行拓扑分析，根据度、最大邻域分量、紧密度、放射性、压力5 种方法筛选出关键基因。

2 结果

2.1 DEGs 的筛选

使用R 软件对GSE55457 数据集中的20 例滑膜样品归一化后，比较各基因在OA 滑膜组织和正常者滑膜组织之间的表达差异，筛选同时符合校正后的P ＜0.05 和|LogFC|>1 的差异基因，共得到203 个DEGs，包括上调基因55 个和下调基因148 个。见图1（封三）。

图1 差异表达分析的火山图

2.2 OA 自噬基因的筛选

通过GeneCards 网站获得符合条件即自噬关联评分≥1 分的自噬基因5786 个，再与DEGs 交集，共获得65 个OA 自噬基因。见图2。

图2 自噬基因与骨关节炎差异基因的交集

2.3 GO 和KEGG 富集分析

利用R 软件对筛选出的65 个OA 自噬基因进行GO 和KEGG 富集分析，结果显示有279 个GO 条目和44 个KEGG 条目；其中与OA 密切相关的GO 分析主要富集在调控DNA 生物合成过程、负调控蛋白质磷酸化、转移酶活性的负调控和血管生成的调节等（图3，见封三）；与OA 高度相关的KEGG 通路有MAPK、TNF、JAK-STAT 和PI3K-Akt 等通路（图4，见封三）。

图3 GO 富集分析结果

图4 KEGG 富集分析结果

2.4 PPI 网络构建与模块分析

通过STRING 数据库构建交集基因的PPI 网络，其节点数目有64 个，边数是129 条，PPI 网络的P 值小于1.0e-16（图5）。另外利用MCODE 插件对网络进行聚类，得到一个较为重要的模块，其包含9 个基因（图6），对其进行KEGG 富集分析得到MAPK、JAKSTAT 和PI3K-Akt 信号通络较为重要（图7，见封三），这与OA 自噬基因富集分析得到结果大致相同。

图5 交集基因的蛋白质-蛋白质相互作用网络

图6 蛋白质-蛋白质相互作用网络中较为重要的模块

图7 模块基因的KEGG 富集分析

2.5 Hub 基因的选择

本研究通过5 种计算方法均筛选出前10 位的关键基因，并利用“Venn”软件对其进行交集和可视化，获得7 个关键基因，包括VEGFA、JUN、CTNNB1、MYC、NR4A1、MAP2K7 和DUSP1。见图8。

图8 5 种计算方法前10 个基因的交集

3 讨论

OA 是最常见的关节退行性疾病，其发病机制尚未明确，而自噬可能在OA 发病中起到重要作用[3]，同时鉴于OA 发病与复杂的细胞因子网络密切相关[6]，因而通过生物信息学挖掘参与OA 发生、发展的自噬重要基因、通路，探究其分子机制。

本研究通过生物信息学分析发现MAPK、JAKSTAT 和PI3K-Akt 通路在OA 自噬中发挥重要作用。在OA 中，MAPK 途径的激活不仅增加软骨和滑液中IL-1、TNF-α、NO 和基质金属蛋白酶-13 的含量[7]，而且可抑制自噬的发生[8]，同时通过抑制MAPK 通路引起下游反应可有效地控制OA 的进展[7]。JAK-STAT 通路的发生是通过激活JAK 激酶而使STAT 磷酸化，接着2 个磷酸化的STAT 在核内与靶基因启动子中特异DNA 序列结合，最终调节细胞增殖、分化、凋亡等[9]；STAT3 作为一种新的自噬调节剂，其通过与PKR 结合使磷酸化的eIF2A 失活以减少自噬途径[10]，因此JAK-STAT 信号通路也介导着自噬的发生[11]，但通过该途径探讨OA 的机制研究甚少，后期可通过实验进一步探讨其相关性。PI3K-AKT 作为调控自噬的经典传导通路，PI3K 的激活可使Akt 磷酸化并引起下游一系列反应，进而导致软骨细胞的炎症反应及自噬相关基因和蛋白表达水平下降[12-13]，而通过抑制PI3K/AKT/mTOR 信号通路可提高细胞自噬来改善软骨的损伤[14]。

此外，本研究通过进一步分析发现了7 个关键基因。VEGFA 可引发滑膜增生和软骨破坏，导致OA 的发生、发展[15-16]。虽然未见有VEGFA 与OA 自噬的报道，但认为由于自噬功能低下而无法抑制VEGFA 表达，可导致肿瘤细胞的持续繁殖，通过激活自噬可有效控制病情发展[17]。JUN 属于原癌基因c-Jun 的表达产物，在转导偶联信息传递中起到核内“第三信使”作用。c-Jun在膝骨关节炎软骨下骨组织中呈高表达状态[18]，可加速OA 软骨损伤。CTNNB1（又称“β-catenin”）是Wnt 信号通路必不可少的蛋白，该途径可调节自噬、细胞凋亡和增殖等[19]。β-catenin 在OA 滑膜组织中呈高表达，不仅可引起滑膜成纤维细胞病理改变[20]，而且可抑制软骨细胞自噬的激活，导致OA 发展[21]。MYC基因家族成员包括c-myc、n-myc 及l-myc，通过影响细胞增殖、分化、自噬和凋亡过程在疾病中发挥作用[22]。c-myc 在OA 软骨中呈高表达，参与OA 发展[23]；有效抑制MYC 表达可增强自噬水平，从而延缓OA 的进程[24]。NR4A1 被认为是调控OA 发展的一种调节剂[25]，其作为TNF/p38 细胞信号传导的下游介质，通过影响软骨细胞的凋亡或自噬在OA 发展中发挥作用[26]。MAP2K7 是MAPK 信号通路中一种重要的激酶，通过JNK 途径可调节OA 中的炎症反应，是研发药物治疗OA 的新靶标[27]。DUSP1（又称“MKP-1”）可通过MAPK通路发挥对OA 的调控作用[28]；另外有研究发现降低或敲除DUSP1 的表达，可促进MAP1/LC3-Ⅱ水平升高，激活细胞自噬[29]，因而DUSP1 可能通过自噬影响OA 进展。

本研究通过对OA 基因表达谱芯片数据与OA自噬基因的分析，挖掘出OA 中与自噬激活相关的重要基因、通路，从自噬角度阐明OA 的发病机制，为研究OA 的分子机制提供一种新的思路。然而，本研究根据大数据预测出的分子机制，仍需进一步通过实验加以验证。