骨关节炎关节滑膜组织差异表达基因生物学功能及与免疫细胞浸润的关系分析

冯重阳，姬振伟，王志学，张智翔，王银歌，2，吴鹏 ，丁勇

1 空军军医大学第二附属医院骨科，西安 710038；2 中国人民解放军联勤保障部队第922医院骨科

骨关节炎（OA）是常见的关节退行性疾病，累及整个关节包括滑膜炎症、软骨下骨硬化、退变和骨赘形成，最终表现为关节疼痛、僵硬和功能障碍［1］。随着全球人口老龄化加剧，OA 患病人数不断增加，影响着约3.8%的人口，造成了巨大的社会经济负担［2］。目前OA 发病机制尚不清楚，多数学者认为，滑膜炎症是OA 的启动环节，是导致患者出现疼痛和临床症状的一个关键因素［3］。AYRAL 等［4］采取连续关节镜检查证实早期OA 存在滑膜炎症，它是OA症状和疾病发展的一个重要病理环节。滑膜炎症由先天免疫细胞参与调控。CD4+和CD8+T 细胞以及先天免疫系统的单核细胞/巨噬细胞的异常都可能导致OA 进行性加重［5-7］，在疾病恶化中起重要作用。因此探究OA 与免疫细胞之间的关系，为揭示OA 发生机制具有重要意义。

生物信息学作为一门新兴学科，被用于快速、大规模地获取生物信息，揭示疾病的分子机制。本研究从Gene Expression Omnibus 数据库（GEO）选取与OA 相关的基因芯片GSE55235、GSE55457、GSE1919和GSE12021，采用Meta 分析的方法筛选OA 和正常样本滑膜组织的差异表达基因（DEGs），对DEGs 进行生物学功能分析，并构建蛋白质—蛋白质相互作用（PPI）网络筛选核心差异表达基因，通过CIBERSORT 算法分析正常和OA 关节滑膜组织之间22 种免疫细胞浸润比例，以及核心差异表达基因与免疫细胞百分比的相关性，为探讨OA 的免疫机制提供理论依据。

1 材料与方法

1.1 OA 数据集及其来源 在 GEO（http：//www.ncbi. nlm. nih. gov/geo）公共数据库中，检索与OA相关的基因芯片，符合以下条件：基因芯片为基因组mRNA 转录组数据的表达谱；包含OA 和正常对照患者滑膜组织的检测数据；排除药物干预或转染的数据集。根据筛选条件获得2 个平台的4个 mRNA 数据集，分别为 GSE55235、GSE55457、GSE1919 和 GSE12021，共有 69 个样本，其中包含35 个 OA 和 34 个 健 康 对 照 样 本 。GSE55235、GSE55457 和 GSE12021 所用平台为 GPL96（Affymetrix Human Genome U133A Array），GSE1919 所 用平 台 为 GPL91（Affymetrix Human Genome U95A Array）。GSE55235 数据集中OA 关节滑膜组织样本 10 个、正常关节滑膜组织 10 个；GSE55457 数据集中OA 关节滑膜组织样本10 个、正常关节滑膜组织10 个；GSE12021 数据集中OA 关节滑膜组织样本 10 个、正常关节滑膜组织 9 个；GSE1919 数据集中OA 关节滑膜组织样本5 个、正常关节滑膜组织5 个。

1.2 OA 与正常关节滑膜组织DEGs 筛选 应用在线分析平台 NetworkAnalyst（https：//www. networkanalyst.ca/）对四个数据集原始数据整合分析［8］。首先将每个数据集上传、注释和处理，而后将四个数据集进行整合归一化处理，采用Fisher's分析方法筛选DEGs（Combined P＜0.05，|logFC|＞2）。对上调和下调前50的DEGs进行层次聚类，使用R软件的pheatmap软件包绘制每个数据集的热图。

1.3 OA 与正常滑膜组织DEGs 的生物学功能分析 GO 分析［9］包括细胞成分（CC）、生物学过程（BP）和分子功能（MF）。KEGG 是一个关于基因组、酶途径和生物化学物质的在线数据库［10］。 应用DAVID 在线生物信息学数据库（https：//david.ncifcrf. gov/，版本6.8）进行基因本体功能及信号通道富集分析，以每个条目里基因数＞10，P＜0.05作为筛选条件。

1.4 OA与正常关节滑膜组织核心差异表达基因筛选 将DEGs 上传至STRING 在线分析平台（https：//string-db. org/，版本 11.5［11］）构建 PPI 网络，使用Cytoscape3.6.8 软件可视化PPI 网络。利用 MCODE 插 件 以 Degree Cut-off =2，Node Score Cutoff =0.2，K-Core =2，Max. Depth = 100 为筛选标准，得到网络中联系最紧密的模块，展示排名前3的模块，选取每个模块中得分最高的基因为核心差异表达基因［12］。

1.5 核心差异表达基因与免疫细胞浸润的相关性分析 通过CIBERSORT算法对NetworkAnalyst分析平台矫正、整合得到的四个数据集的芯片数据进行计算，得到22 种免疫细胞在每个样本中所占比例。使用vioplot 软件包绘制小提琴图。使用ggstatsplot软件包对筛选的核心差异表达基因水平与22 种免疫细胞百分比进行Spearman相关性分析。

2 结果

2.1 OA 与正常关节滑膜组织的DEGs 及其生物学功能 共筛选得到1 418 个OA 与正常滑膜组织的DEGs，其中上调 692 个，下调 726 个。GO 功能分析结果显示，DEGs 主要富集在酰胺结合、DNA 结合转录因子结合、糖胺聚糖结合、整合素结合、核受体结合等。KEGG 分析结果显示，DEGs 主要参与人类T细胞白血病病毒感染、TNF信号通路、类风湿性关节炎、癌症中的蛋白聚糖、爱泼斯坦-巴尔病毒感染等。

2.2 OA 与正常关节滑膜组织核心差异表达基因通过Cytoscape 软件构建了一个包含820 个节点、2 212 条边的 PPI 网络图。RPL37A、MYC 和细胞分裂周期蛋白20（CDC20）可能是OA 的核心差异表达基因。

2.3 OA和正常关节滑膜组织核心差异表达基因与免疫细胞浸润的相关性 与正常关节滑膜组织相比，OA 患者关节滑膜组织中记忆B细胞、浆细胞、辅助性T 细胞和静息肥大细胞所占比例较多，而静息记忆性CD4+T 细胞、嗜酸性粒细胞和活化肥大细胞所占比例较少（P 均＜0.05）。核心差异表达基因与免疫细胞的相关性分析显示，RPL37A 与γδT 细胞、M2 巨噬细胞、单核细胞百分比呈正相关（r 分别为0.36、0.29、0.25，P 均＜0.05），与初始 B 细胞、调节性 T 细胞、CD8+T 细胞、记忆 B 细胞、活化肥大细胞、静息树突状细胞百分比呈负相关（r 分别为-0.38、-0.40、-0.37、-0.32、-0.35、-0.25，P 均＜0.05）。MYC 与静息记忆性T 细胞、嗜酸性粒细胞、活化树突状细胞、单核细胞、活化肥大细胞百分比呈正相关（r分别为0.52、0.37、0.35、0.34、0.29，P均＜0.05），与记忆 B 细胞、调节性 T 细胞、CD8+T 细胞、静息树突状细胞、辅助性T 细胞百分比呈负相关（r 分别为-0.49、-0.48、-0.4、-0.27、-0.27，P 均＜0.05）。CDC20 与 M2 巨噬细胞、静息肥大细胞、辅助性T 细胞百分比呈正相关（r 分别为0.59、0.30、0.25，P 均＜0.05），与活化肥大细胞、嗜酸性粒细胞、初始B 细胞、活化树突状细胞、静息记忆性T 细胞百分比呈负相关（r 分别为-0.40、-0.35、-0.34、-0.26、-0.25，P均＜0.05）。

3 讨论

OA 是常见的关节疾病，主要累及全身大关节，其中膝关节是最常见的累及部位，影响患者的日常生活，导致残疾［13］。但目前还没有可以有效预防或逆转关节损伤的治疗手段，终末期患者只能行关节置换。近年，越来越多的证据表明，OA 的早期到晚期始终伴随着滑膜炎症［14］，且滑膜炎症会对邻近的骨和软骨造成损害，加重疾病进展［4，14］。此外，滑膜炎症中免疫细胞浸润在OA 的发生发展中也起着重要作用［15］。因此，探讨OA 与机体免疫之间的关系，有助于进一步明确OA发病机制。

本研究通过筛选得到四个基因芯片，使用NetworkAnalyst 分析平台整理、合并四个OA 关节滑膜组织基因表达芯片，随后得到1 418个DEGs，其中上调692个、下调726个。GO富集分析结果显示DEGs主要涉及酰胺结合、DNA 结合转录因子结合、糖胺聚糖结合、整合素结合、核受体结合等。KEGG 富集通路显示主要涉及人类T 细胞白血病病毒感染、TNF 信号通路、类风湿性关节炎、癌症中的蛋白聚糖、爱泼斯坦-巴尔病毒感染等相关通路。GO 富集分析表明，差异基因的功能与OA 的调控密切相关，糖胺聚糖的减低会导致软骨细胞的减少以及组织的降解，致使关节功能丧失［16］，整合素作为软骨细胞表面受体，介导软骨细胞与细胞外基质的黏附，多种整合素及其配体是OA 疾病进展的生物标志物［17］。KEGG 富集分析表明，差异基因的相关通路除了与病毒感染性免疫反应相关外，在类风湿性关节炎、TNF 信号通路中基因富集较多，这些通路参与OA滑膜炎症因子渗出、免疫细胞浸润的过程［18］。TNF信号通路在OA 的发生机制中已被广泛研究，它可通过调控NF-κB 和MAPK 两条关键信号通路减少炎症因子的分泌，从而延缓OA的进展［19］。

本研究通过构建DEGs 的PPI 网络，利用MCODE 插件对PPI 进行分析，根据得分筛选出RPL37A、MYC、CDC20 为 OA 发展过程中的核心差异表达基因。RPL37A 是编码核糖体蛋白的基因［20］。本研究显示RPL37A 在正常人中高表达，与MI等［21］的分析结果一致。女性一直被认为是OA 发生的危险因素，YANG 等［22］的研究发现在女OA患者中 RPL37A 表达降低，确定了 RPL37A 是女性 OA 的核心基因。在斑马雀中，雌性表达RPL37A 的细胞总数低于雄性，并且随着年龄的增长，雌性体内RPL37A 反而下降，这一结果有可能解释为什么OA经常发生在老年女性［23］。MYC（也称为c-Myc）属于MYC 原癌基因家族，是一种转录因子，定位于8q24，用来调控细胞周期、生长、凋亡、分化、代谢等核心过程［24］。FISCH 等［25］研究确定了 MYC 可能是 OA 发病的核心基因。MYC 通过消除人体内和体外肾器官中的炎症和纤维化，从而降低炎症反应［26］。MAPK信号通路参与了人类OA 软骨细胞的凋亡和促炎症细胞因子的表达［27］，其与MYC 密切相关。WU 等［28］研究证明，MYC 可以通过抑制MAPK 信号通路，抑制软骨细胞的凋亡，促进软骨细胞的增殖，来阻止OA 的进展，本研究与文献中结果一致，证实了MYC在OA 患者滑膜组织样本中表达降低，致使炎症发生。CDC20 通常被认为是致癌因子，主要通过活化APC 形成一种E3 泛素连接酶复合物APCCDC20，参与其下游底物的降解过程，来调节细胞周期、促进细胞凋亡等，与肿瘤的发生关系密切。本研究发现，CDC20 在OA 关节滑膜组织中高表达，提示CDC20可能在OA 发生中起重要作用，这在以往研究中未见报道。其原因可能是CDC20 与细胞周期密切相关，而细胞周期的异常调节与OA 的发展密切相关［29］。

慢性滑膜炎的免疫细胞浸润以巨噬细胞为标志。OA 患者滑膜中T 细胞和B 细胞的含量高于健康患者［3］。ROSSHIRT 等［15］研究表明，OA 出现免疫细胞浸润，包括CD4+T 细胞、CD8+T 细胞。为了进一步了解OA 滑膜炎症中免疫细胞浸润情况，我们使用CIBERSORT 算法进行了全面评估，CIBERSORT 是一个利用RNA-seq 数据评估免疫细胞表达的分析工具，可以用于深入了解正常和OA 关节滑膜组织中的免疫微环境。通过对比发现与正常关节滑膜组织相比，OA 患者关节滑膜组织中记忆B 细胞、浆细胞、辅助性T 细胞和静息肥大细胞比例较高，而静息记忆性CD4+T 细胞、嗜酸性粒细胞和活化肥大细胞比例较少，与GE 等［30］部分结果一致。当滑膜组织中静息肥大细胞的浸润比例增高时，可导致OA 患者的软骨结构发生损伤［31］。辅助性T 细胞会产生分解代谢细胞因子包括IL-2、IFN-γ、TNF-α，随着滑膜组织中浸润过多的辅助性T 细胞，其分泌的分解代谢细胞因子也会随之增加，导致OA的进展［32］。DE LANGE-BROKAAR 等［3］发现，浆细胞也可能影响着OA的进展。

我们进一步分析了核心差异表达基因与免疫细胞浸润的关系，发现RPL37A与γδT细胞、M2巨噬细胞、单核细胞呈正相关，与初始B 细胞、调节性T 细胞、CD8+T细胞、记忆B细胞、活化肥大细胞、静息树突状细胞呈负相关。MYC 与静息记忆性CD4+T 细胞、嗜酸性粒细胞、活化树突状细胞、单核细胞、活化肥大细胞所占比例呈正相关，与记忆B 细胞、调节性T 细胞、CD8+T 细胞、静息树突状细胞、滤泡辅助性T 细胞所占比例呈负相关。CDC20 与M2 巨噬细胞、静息肥大细胞、滤泡辅助性T 细胞呈正相关，活化肥大细胞、嗜酸性粒细胞、初始B 细胞、活化树突状细胞、静息记忆性CD4+T 细胞所占比例呈负相关。说明核心差异表达基因与免疫细胞的浸润关系密切。

综上所述，OA 关节滑膜组织DEGs 共1 418个，其中692 个上调、726 个下调，基因本体功能主要有酰胺结合、DNA 结合转录因子结合等，参与人类T 细胞白血病病毒感染、TNF 信号通路等；OA 的核心差异表达基因是RPL37A、MYC、CDC20；记忆B 细胞、浆细胞等可能参与OA 关节滑膜炎症的调控，RPL37A、MYC、CDC20 与免疫细胞浸润关系密切。本研究加深了对OA 关节滑膜炎症免疫分子机制的认识，可能为OA 早期诊断及治疗提供潜在标志物，然而，这些还需要临床以及进一步实验验证。