楊延輯，廖文波

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

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

1 材料與方法

1.1 數(shù)據(jù)的下載及處理 本研究中從公共數(shù)據(jù)庫GEO(GENE EXPRESSION OMNIBUS,GEO)下載骨關(guān)節(jié)炎滑膜相關(guān)芯片數(shù)據(jù)集。篩選條件為：①骨關(guān)節(jié)炎；②人類；③滑膜組織；④無藥物及手術(shù)干預(yù)?；诤Y選條件下載表達(dá)數(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ù)芯片進(jìn)行預(yù)處理，將GSE1919，GSE55235兩個數(shù)據(jù)集作為訓(xùn)練集進(jìn)行數(shù)據(jù)合并，將GSE82017,GSE55457作為兩個獨立驗證集，對樣本中都不表達(dá)的探針、存在的缺失值或基因與多個探針存在對應(yīng)關(guān)系等特殊情況，進(jìn)行標(biāo)準(zhǔn)化、探針過濾、缺失值填充以及探針合并等，并對合并數(shù)據(jù)表達(dá)矩陣進(jìn)行l(wèi)og2對數(shù)轉(zhuǎn)換，對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理。

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

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

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

2 結(jié)果

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

(紅色高表達(dá)，藍(lán)色低表達(dá))

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

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

2.3 機器學(xué)習(xí)篩選OA標(biāo)志物 對獲得的474差異基因進(jìn)一步通過機器學(xué)習(xí)的方式篩選特征基因，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個特征標(biāo)志基因BCL6、DDIT4、SCRG1，繪制韋恩圖(見圖4)。

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

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

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

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