楊延輯，廖文波

(1.遵義醫(yī)科大學附屬醫(yī)院 脊柱外科，貴州 遵義 563099；2.遵義醫(yī)科大學第二附屬醫(yī)院 骨外科，貴州 遵義 563000)

骨性關(guān)節(jié)炎(Osteoarthritis,OA)是由多種因素引起的以關(guān)節(jié)軟骨退行性變,當骨關(guān)節(jié)炎發(fā)展到一定程度時形成不可逆疾病，世界上超過 10% 的人口患有OA，65歲以上人群患病率更高，OA 造成的經(jīng)濟和社會負擔正在迅速增加，并嚴重影響患者的生活質(zhì)量[1]。關(guān)節(jié)置換是關(guān)節(jié)炎目前終末期的主要治療手段，其中膝骨關(guān)節(jié)炎(Knee osteoarthritis,KOA)患病率高、花費高是導致世界殘疾的第四大原因[2]，目前臨床診斷大多靠影像學、檢驗學、查體等；隨著社會人口的老齡化 ,該病的發(fā)生率越來越高，目前的藥物治療大多是對癥的，迄今為止還缺乏改善疾病的 OA 藥物[3]；因此骨關(guān)節(jié)炎早期診斷、及時防治有其重要現(xiàn)實意義[4]。骨關(guān)節(jié)炎的發(fā)病機制目前尚不明確，以往骨關(guān)節(jié)炎的研究熱點大多集中于軟骨組織，現(xiàn)研究發(fā)現(xiàn)滑膜、成纖維樣滑膜細胞和軟骨下骨細胞也參與了骨關(guān)節(jié)炎的發(fā)病過程[5]。本研究擬通過生物信息學方法，對GEO公共數(shù)據(jù)庫中有關(guān)OA滑膜的基因表達譜芯片數(shù)據(jù)進行差異表達基因分析、功能富集、嘗試用機器學習的方式篩選骨關(guān)節(jié)炎關(guān)鍵基因，探討OA發(fā)生發(fā)展中調(diào)節(jié)基因，為探索其診斷和治療靶點提供生物學信息依據(jù)。

1 材料與方法

1.1 數(shù)據(jù)的下載及處理 本研究中從公共數(shù)據(jù)庫GEO(GENE EXPRESSION OMNIBUS,GEO)下載骨關(guān)節(jié)炎滑膜相關(guān)芯片數(shù)據(jù)集。篩選條件為：①骨關(guān)節(jié)炎；②人類；③滑膜組織；④無藥物及手術(shù)干預。基于篩選條件下載表達數(shù)據(jù)集GSE1919，GSE55235，GSE82017，GSE55457；GSE1919包含5個OA 滑膜組織樣本和5個正?；そM織樣本，GSE55235包含10個OA滑膜組織樣本和10個正?；そM織樣本，GSE82017包含10個OA滑膜組織樣本和7個正?；そM織樣本，GSE55457包含10個OA滑膜組織樣本和10個正?；そM織樣本，利用R軟件讀取下載相關(guān)基因的原始數(shù)據(jù)，對數(shù)據(jù)芯片進行預處理，將GSE1919，GSE55235兩個數(shù)據(jù)集作為訓練集進行數(shù)據(jù)合并，將GSE82017,GSE55457作為兩個獨立驗證集，對樣本中都不表達的探針、存在的缺失值或基因與多個探針存在對應關(guān)系等特殊情況，進行標準化、探針過濾、缺失值填充以及探針合并等，并對合并數(shù)據(jù)表達矩陣進行l(wèi)og2對數(shù)轉(zhuǎn)換，對數(shù)據(jù)進行標準化處理。

1.2 差異基因及功能富集 利用R軟件(4.1.2版本)Limma包進行差異表達分析，設定P<0.05,基因表達差異倍數(shù)(FoId change，F(xiàn)c)絕對值≥1為篩選條件，篩選獲得差異表達基因(Differentially expressed gene,DEG)。利用cluster-Profiler、org.Hs.eg.db及enrichplot包對獲得的差異基因進行基因本體論(Gene Ontology ,GO)和京都基因與基因組百科全書(Kyoto encyclopedia of genes and genomes,KEGG)富集分析，以P<0.05為篩選標準；并使用ggplot2包進行可視化處理。

1.3 機器學習篩選特征基因 機器學習(Machine learning,ML)是一門研究計算機如何模擬人類進行學習數(shù)據(jù)分析的科學[6]，在處理高維度、大批量數(shù)據(jù)等方面較傳統(tǒng)方法有顯著優(yōu)勢，目前已在臨床獲得廣泛應用，對上述獲得的差異基因通過R軟件使用lasso回歸算法[7]和SVM-RFE 支持向量機遞歸特征消除算法[8]兩種機器學習方式篩選特征基因，LASSO 是一種回歸算法，SVM-RFE 由Guyon等[9]提出，用于癌癥分類中基因的選擇，在臨床已得到應用，如自閉癥識別[10]、前列腺組織病理學分級等[11]。對兩種方式獲得的基因結(jié)果取交集，應用R軟件“venneuler”包繪制韋恩圖對結(jié)果進行可視化。

1.4 診斷標志物的驗證 對上述取交集獲得的特征基因，在GSE82017,GSE55457兩個獨立的數(shù)據(jù)集中來驗證獲得的基因作為診斷標志物的價值，通過繪制接受者操作特征曲線(Receiver operating characteristic，ROC)曲線評價其診斷價值，以P<0.05為閾值來確定。計算ROC曲線下面積(Area under the ROC，AUC)，AUC的取值在0～1，AUC越大，說明預測性能越好。

2 結(jié)果

2.1 差異基因篩選 對GSE1919及GSE55235兩個數(shù)據(jù)進行合并，矯正、標準化處理后得到8 920個基因，進行差異分析共獲得474個差異基因，其中226個表達上調(diào)基因，248個表達下調(diào)基因，繪制火山圖并對表達排名前100的差異基因繪制熱圖(見圖1)。

(紅色高表達，藍色低表達)

2.2 GO及KEGG富集結(jié)果 GO富集分析發(fā)現(xiàn)生物過程(Biological process,BP)主要富集在正向調(diào)節(jié)白細胞的激活，正向調(diào)節(jié)細胞激活、附著力，細胞趨化性，對脂多糖的反應，淋巴細胞活化的正調(diào)節(jié)等，分子功能(Molecular function,MF)中含膠原蛋白細胞外基質(zhì)、等離子外側(cè)膜、膜筏、膜微區(qū)、胞吞泡富集明顯，細胞組分(Cellular component,CC)顯示與受體配體活性、信號受體激活劑活性、糖胺聚糖結(jié)合、硫化合物結(jié)合等相關(guān)；KEGG富集分析顯示其在脂質(zhì)和動脈粥樣硬化、1型人類 T 細胞白血病病毒感染、類風濕關(guān)節(jié)炎、MAPK信號通路、細胞-細胞因子受體相互作用、NF-kappa B 信號通路、破骨細胞分化、AGE-RAGE 信號通路、腫瘤壞死因子信號通路、IL-17信號通路、趨化因子信號通路等(見圖2)。

圖2 GO、KEGG富集分析結(jié)果

2.3 機器學習篩選OA標志物 對獲得的474差異基因進一步通過機器學習的方式篩選特征基因，lasso算法選擇lambda.min參數(shù)[12]獲得14個特征基因分別為BCL6、DDIT4、KLF9、GADD45A、SIK1、HNRNPA1、MRC2、MTHFD2、KDELR3、CX3CR1、SCRG1、SLC2A3、ABCC3、ABL2，通過SVM-RFE算法得到19個特征基因分別為SELL、SCRG1、HLA-DMB、BCL6、UCP2、HLA-DMA、MIR6883、MXRA5、DDIT4、NFIL3、TNFAIP3、SORL1、MIR8071-2、STC1、NEDD9、KLF4、FAM107A、SNORD10、MAFF(見圖3)，將兩種方式獲得的基因取交集，獲得3個特征標志基因BCL6、DDIT4、SCRG1，繪制韋恩圖(見圖4)。

圖3 機器學習Lasso結(jié)果及VSM-RFE結(jié)果

圖4 兩種機器學習結(jié)果取交集得到的3個基因

2.4 特征基因的驗證價值 對獲得的BCL6，DDIT4，SCRG1特征基因，在GSE82017,GSE55457兩個獨立的數(shù)據(jù)集中進行外部驗證3個基因的診斷價值。ROC曲線結(jié)果顯示BCL6在兩個數(shù)據(jù)集的AUC值分別為0.9及0.83；SCRG1在兩個數(shù)據(jù)集的AUC值分別為0.829及0.93；DDIT4的AUC值分別為0.757及0.7，3個診斷基因AUC值均>0.5，結(jié)果顯示其均有較高診斷價值(見圖5)。

A:GSE82107中驗證結(jié)果;B:GSE55457中驗證結(jié)果。

3 討論

骨性關(guān)節(jié)炎是退行性關(guān)節(jié)疾病，多因素復雜作用所致；目前仍缺乏一種能夠得到廣泛認可的、臨床上普遍使用的OA診斷的生物標記物。為了探索骨關(guān)節(jié)炎的診斷標志物，本研究通過從公共數(shù)據(jù)庫GEO下載數(shù)據(jù)集GSE1919，GSE55235，GSE82017，GSE55457；共包含35個OA滑膜組織樣本和，32個正常滑膜組織，將GSE1919，GSE55235作為訓練集合并數(shù)據(jù)后篩選獲得DEGs；對DEGs進行富集分析，其富集結(jié)果中現(xiàn)已有大量研究證實核轉(zhuǎn)錄因子kB(Nuclear factor-kapa B,NF-kB)信號通路的激活在OA的發(fā)生，發(fā)展中起著重要作用; NF-kB信號通路的激活增加炎癥因子的釋放，導致關(guān)節(jié)軟骨的破壞，以及滑膜炎的發(fā)生[13]，MAPK 信號通路在早期 OA 中的調(diào)節(jié)作用，抑制MAPK信號通路可緩解OA的進展[14]，使用lasso回歸算法和SVM-RFE 支持向量機遞歸特征消除算法得到的結(jié)果取交集獲得3個特征基因，分別為BCL6，DDIT4，SCRG1。BCl6是位于3q27染色體的轉(zhuǎn)錄抑制因子[15]，其參與細胞的增殖、分化、凋亡及炎癥過程等[16]。目前研究發(fā)現(xiàn)，BCL6 表達水平的升高提示患者滑膜炎癥加重狀態(tài)，在類風濕關(guān)節(jié)炎患者中BCL6的變化與患者病情程度及膜炎病理評分呈正相關(guān)[17]。90%的骨關(guān)節(jié)炎患者存在滑膜病變，且滑膜的病變程度與關(guān)節(jié)嚴重的疼痛和功能障礙有關(guān)[18]，滑膜炎癥可促進促炎因子和疼痛神經(jīng)遞質(zhì)的產(chǎn)生[19]。此外，滑膜炎的發(fā)生還可能促進軟骨退變，滑膜炎通過釋放促炎介質(zhì)和軟骨破壞因子引起軟骨損傷，從而在OA中發(fā)揮引發(fā)作用[20]；Hou等[21]研究證實，BCL6高表達可以抑制NFATc1的活化來調(diào)節(jié)RANKL，而在OA中OPG/RANK/RANKL 調(diào)節(jié)系統(tǒng)的功能障礙與軟骨下骨的組織學改變之間存在密切相關(guān)性，RANKL 不僅參與影響軟骨下骨的通路，而且在骨關(guān)節(jié)炎軟骨中表達，RANKL 的產(chǎn)生會導致軟骨破壞[22]，并且XU等[23]通過實驗證實，抑制RANKL 誘導的破骨細胞生成，能改善軟骨細胞炎癥減緩OA的進展，綜上推測BCL6可能在OA起著一定作用，或可成為潛在治療靶點，但目前BCL6大多集中在腫瘤，在其他疾病領(lǐng)域的研究近期才逐漸增多，仍需進一步臨床實驗研究來證實。目前已有研究發(fā)現(xiàn)，SCRG1在人關(guān)節(jié)軟骨中特異性表達[24]，間充質(zhì)干細胞(MSCs)是非造血基質(zhì)細胞，具有自我更新和分化成間充質(zhì)細胞的能力[25]，研究發(fā)現(xiàn)SCRG1正向調(diào)節(jié)hMSC自我更新、遷移和成骨分化和成軟骨分化[26]，SCRG1是一種促進軟骨基因表達的轉(zhuǎn)錄因子，ADAMTS9-AS2高表達可以逆轉(zhuǎn) miR-942-5p對SCRG1的抑制[27]。新近研究發(fā)現(xiàn)，SCRG1通過 Wnt5a 促進臍帶間充質(zhì)干細胞的成軟骨分化[28]；骨關(guān)節(jié)炎可導致軟骨撕裂和軟骨細胞丟失，軟骨變薄，磨損退化[4]，而SCRG1在人關(guān)節(jié)軟骨中特異性表達，其可能在OA的發(fā)生發(fā)展中起到一定作用，但目前整體關(guān)于SCRG1在各個疾病領(lǐng)域的研究還相對較少，其在骨關(guān)節(jié)炎中的具體作用機制還需進一步研究。DDIT4是一種腫瘤相關(guān)蛋白，也是一種代謝和免疫相關(guān)蛋白，應激反應蛋白，可通過mTOR 途徑調(diào)節(jié)蛋白質(zhì)合成、自噬和細胞凋亡；在化療、缺氧和DNA損傷等應激條件下高度表達[29-30],DDIT4也被稱為發(fā)育和DNA損傷反應調(diào)控(REDD1)，關(guān)節(jié)軟骨退化是OA發(fā)生的主要因素之一，涉及mTOR信號傳導增加和自噬缺陷，REDD1 在正常人關(guān)節(jié)軟骨中高度表達，并在衰老和 OA 期間減少，REDD1是mTOR的抑制劑；通過建立OA實驗小鼠研究發(fā)現(xiàn)，REDD1缺乏會增加小鼠實驗性OA的嚴重程度，并且在滑膜、半月板和軟骨下骨中觀察到類似的差異[31]，另有研究發(fā)現(xiàn)REDD1 在正常人關(guān)節(jié)軟骨中高度表達，并在衰老和 OA 期間減少，通過實驗觀察發(fā)現(xiàn)REDD1在人和小鼠軟骨以及小鼠半月板和滑膜中的表達隨著年齡的增長而降低[32]。通過從鐵死亡相關(guān)基因從FerrDb數(shù)據(jù)庫(http://zhounan.org/ferrdb)[33]查詢DDI4是鐵死亡相關(guān)基因；研究發(fā)現(xiàn)REDD1 表達被發(fā)現(xiàn)是許多病理的早期生物標志物，包括炎癥性疾病，REDD1 的表達與細胞凋亡、活性氧(ROS)的堆積和導致組織損傷的炎癥激活有關(guān)，其可通過NF-κB、一氧化氮合酶和超氧化物歧化酶/谷胱甘肽過氧化物酶/NAPDH 氧化酶途徑激活炎癥和 ROS 產(chǎn)生[34]。并且新進研究發(fā)現(xiàn)，mTORC1作為mTOR的復合物，實驗證明鐵螯合通過多種途徑抑制 mTORC1，而鐵對 mTORC1 的激活至關(guān)重要，其中證實REDD1 通路部分參與鐵螯合誘導的 mTORC1 抑制[35]；鐵死亡作為是近年發(fā)現(xiàn)的一種新的程序性細胞死亡類型，這種新型死亡方式與活性氧(ROS)及細胞內(nèi)鐵相關(guān)，鐵螯合劑可以抑制這一過程[36]，現(xiàn)在越來越多的研究發(fā)現(xiàn)鐵死亡與OA的發(fā)生、發(fā)展相關(guān)，并且認為鐵代謝異常是OA發(fā)展的影響因素之一[37]，綜上,DDIT4可能成為臨床在鐵死亡方向作為OA治療及診斷的一個新標志物。

綜上所述，在本研究中使用綜合生物信息學分析和機器學習方法篩選出BCL6，SCRG1，DDIT4鑒定為OA相關(guān)的特征基因，它們可能在OA的發(fā)生發(fā)展中充當重要的調(diào)節(jié)因子；并且在2個獨立驗證集ROC曲線中AUC值來驗證3個相關(guān)基因，均獲得較好的診斷價值，其中BCL6，SCRG1較DDIT4具有更高的診斷價值，DDIT4作為鐵死亡相關(guān)基因，本研究通過生物信息學分析也證明鐵死亡相關(guān)基因在OA發(fā)生發(fā)展中的作用，為OA的發(fā)病機制和治療提供新的見解，但本研究是基于生物信息學及機器學習方法進行分析與解讀，還需進一步的實驗來進一步驗證。