基于生物信息學方法鑒別骨關節炎自噬的關鍵基因和途徑

袁長深 容偉明 李 哲 官巖兵 梅其杰 段 戡

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 的分子機制提供一種新的思路。然而，本研究根據大數據預測出的分子機制，仍需進一步通過實驗加以驗證。