骨關節炎關節滑膜組織差異表達基因生物學功能及與免疫細胞浸潤的關系分析

馮重陽，姬振偉，王志學，張智翔，王銀歌，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 早期診斷及治療提供潛在標志物，然而，這些還需要臨床以及進一步實驗驗證。