龍通華，張紅參，許冬梅，曾志華，傅子原，關(guān)皓鄴

(1.右江民族醫(yī)學(xué)院研究生學(xué)院，2.右江民族醫(yī)學(xué)院附屬醫(yī)院康復(fù)醫(yī)學(xué)科，3.廣西高校重點(diǎn)實(shí)驗(yàn)室右江流域特色民族藥研究重點(diǎn)實(shí)驗(yàn)室，廣西百色 533000)

骨關(guān)節(jié)炎(osteoarthritis，OA)是最常見(jiàn)的慢性退行性骨骼系統(tǒng)疾病之一，多見(jiàn)于中老年人。骨關(guān)節(jié)炎的病理特征是軟骨退化、滑膜發(fā)炎、軟骨下改變及骨贅形成，最終導(dǎo)致關(guān)節(jié)功能喪失[1]，因此骨關(guān)節(jié)炎是導(dǎo)致老年人關(guān)節(jié)疼痛、身體殘疾和生活質(zhì)量下降的主要原因之一[2]。目前，骨關(guān)節(jié)炎的發(fā)病機(jī)制仍不清楚，缺少可用于診斷的生物標(biāo)志物，因此，利用生物信息學(xué)方法探索骨關(guān)節(jié)炎的生物標(biāo)志物及其潛在的作用機(jī)制對(duì)于骨關(guān)節(jié)炎的診斷、治療和預(yù)后至關(guān)重要。

有研究發(fā)現(xiàn)，骨關(guān)節(jié)炎的發(fā)病機(jī)制與基因表達(dá)調(diào)控有關(guān)[3-4]?；蛐酒夹g(shù)、基因測(cè)序技術(shù)和生物信息學(xué)分析在基因組水平的數(shù)據(jù)分析中被廣泛應(yīng)用，可識(shí)別與骨關(guān)節(jié)炎發(fā)生發(fā)展相關(guān)的差異表達(dá)基因(differentially expressed genes，DEGs)和功能通路。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis，WGCNA)是一種基于高通量基因表達(dá)譜的算法[5]。它被廣泛用于識(shí)別各種疾病的基因共表達(dá)網(wǎng)絡(luò)，以揭示基因之間的相關(guān)性，并找到疾病顯著相關(guān)的基因模塊[6-8]。基因表達(dá)綜合數(shù)據(jù)庫(kù)(gene expression omnibus，GEO, http: //www.ncbi.nlm.nih.gov/geo)是一個(gè)高通量微陣列和下一代測(cè)序序列功能基因組數(shù)據(jù)庫(kù)[9]?；虮倔w論(gene ontology，GO)是注釋基因和分析基因生物學(xué)過(guò)程的主要生物信息學(xué)工具[10]。京都基因與基因組百科全書(Kyoto encyclopedia of genes and genomes，KEGG)是一個(gè)整合了基因組、化學(xué)和系統(tǒng)功能信息的數(shù)據(jù)庫(kù)，并對(duì)基因的功能進(jìn)行系統(tǒng)化分析[11]。

本研究基于生物信息學(xué)分析方法，通過(guò)GEO數(shù)據(jù)庫(kù)下載含有20個(gè)骨關(guān)節(jié)炎和18個(gè)正常對(duì)照樣本的高通量測(cè)序數(shù)據(jù)集GSE114007進(jìn)行差異分析，利用WGCNA將差異基因分為不同的基因模塊，并篩選出與骨關(guān)節(jié)炎顯著相關(guān)的基因模塊進(jìn)行GO和KEGG富集分析，尋找其潛在的作用機(jī)制。蛋白互作(protein-protein interaction，PPI)網(wǎng)絡(luò)篩選出DDIT3和BCL6作為關(guān)鍵基因。最后從GEO數(shù)據(jù)庫(kù)下載GSE57218和GSE169077表達(dá)譜數(shù)據(jù)對(duì)訓(xùn)練集GSE114007分析結(jié)果進(jìn)行驗(yàn)證。

1 材料與方法

1.1 數(shù)據(jù)下載和預(yù)處理本研究從GEO數(shù)據(jù)庫(kù)搜索并下載與骨關(guān)節(jié)炎相關(guān)的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)和表達(dá)譜數(shù)據(jù)，得到GSE114007、GSE57218和GSE169077三個(gè)數(shù)據(jù)集，其中GSE114007作為訓(xùn)練集，GSE57218和GSE169077作為驗(yàn)證集。GSE114007是轉(zhuǎn)錄組測(cè)序的counts數(shù)據(jù)，有20個(gè)骨關(guān)節(jié)炎樣本和18個(gè)正常對(duì)照樣本，GSE57218有33個(gè)骨關(guān)節(jié)炎樣本和7個(gè)正常對(duì)照樣本， GSE169077有6個(gè)骨關(guān)節(jié)炎樣本和5個(gè)正常對(duì)照樣本。利用R語(yǔ)言軟件對(duì)表達(dá)譜數(shù)據(jù)進(jìn)行l(wèi)og2轉(zhuǎn)換和歸一化處理，根據(jù)平臺(tái)注釋信息將探針I(yè)D轉(zhuǎn)換基因名，數(shù)據(jù)中多個(gè)基因探針I(yè)D對(duì)應(yīng)單個(gè)基因名，將數(shù)據(jù)進(jìn)行求和取最大值處理。

1.2 篩選差異基因使用R語(yǔ)言軟件的“DESeq2”包[12-13]對(duì)GSE114007數(shù)據(jù)的骨關(guān)節(jié)炎樣本與正常對(duì)照樣本進(jìn)行差異分析，篩選DEGs。DEGs篩選標(biāo)準(zhǔn)為矯正P(adj.P.Val)<0.05和|log2FC|>1。

1.3 構(gòu)建DEGs共表達(dá)基因模塊利用“WGCNA”包[5]對(duì)DEGs構(gòu)建共表達(dá)網(wǎng)絡(luò)，首先pickSoftThreshold函數(shù)計(jì)算出相關(guān)系數(shù)為0.86時(shí)對(duì)應(yīng)的軟閾值為10。以軟閾值為10將差異表達(dá)基因構(gòu)建成鄰接矩陣，之后轉(zhuǎn)變?yōu)橥負(fù)渲丿B矩陣[14]，以識(shí)別基因模塊?；蚰K的最小基因數(shù)為30，合并相似基因模塊的切割高度的閾值設(shè)為0.25[5]。

1.4 GO和KEGG富集分析為了進(jìn)一步分析與OA最相關(guān)的基因模塊的基因功能，我們使用“clusterProfiler”包[15]分別對(duì)與OA最相關(guān)的基因模塊的基因進(jìn)行GO富集分析和KEGG通路富集分析，篩選條件為P<0.05。

1.5 PPI分析和關(guān)鍵基因的篩選通過(guò)STRING在線數(shù)據(jù)11.5(https://string-db.org/cgi/input.pl)構(gòu)建與骨關(guān)節(jié)炎最相關(guān)的基因模塊基因的蛋白質(zhì)-蛋白質(zhì)相互作用網(wǎng)絡(luò)關(guān)系，以識(shí)別蛋白質(zhì)之間的相互作用，最低互動(dòng)分?jǐn)?shù)為0.9。并通過(guò)Cytoscape軟件v3.8.2版本進(jìn)行可視化PPI網(wǎng)絡(luò)[16]，利用cytoHubba插件的Degree算法篩選出前十個(gè)基因作為關(guān)鍵基因[17]。

1.6 關(guān)鍵基因的驗(yàn)證在10個(gè)關(guān)鍵基因中有兩個(gè)關(guān)鍵基因(DDIT3和BCL6)在骨關(guān)節(jié)炎中尚未見(jiàn)相關(guān)研究報(bào)道。利用t檢驗(yàn)計(jì)算出關(guān)鍵基因在GSE57218和GSE169077表達(dá)譜數(shù)據(jù)中OA組和正常組的表達(dá)水平，驗(yàn)證關(guān)鍵基因(DDIT3和BCL6)表達(dá)水平的差異。

1.7 統(tǒng)計(jì)學(xué)方法根據(jù)數(shù)據(jù)分布特點(diǎn)，采用非參數(shù)檢驗(yàn)或t檢驗(yàn)分析兩組間差異的統(tǒng)計(jì)學(xué)意義。所有分析均使用R軟件4.0.5版本進(jìn)行分析，adj.P.Val或P<0.05為有統(tǒng)計(jì)學(xué)意義。

2 結(jié) 果

2.1 差異表達(dá)基因篩選在骨關(guān)節(jié)炎和正常對(duì)照兩組中共篩選出1752個(gè)DEGs，其中上調(diào)基因927個(gè)，下調(diào)基因825個(gè)，DEGs的結(jié)果展示如圖1。

2.2 DEGs共表達(dá)網(wǎng)絡(luò)分析我們選取相關(guān)系數(shù)為0.86時(shí)對(duì)應(yīng)的軟閾值為10(如圖2A所示)構(gòu)建無(wú)尺度網(wǎng)絡(luò)進(jìn)行WGCNA分析。通過(guò)WGCNA分析將DEGs構(gòu)建了6個(gè)不同的基因模塊，分別是棕色基因模塊、藍(lán)色基因模塊、藍(lán)綠色基因模塊、綠色基因模塊、黃色基因模塊及灰色基因模塊，如圖2B所示。根據(jù)圖3A所示的基因模塊與骨關(guān)節(jié)炎相關(guān)性分析表明，這些基因模塊都與骨關(guān)節(jié)炎的發(fā)生發(fā)展相關(guān)。圖3B顯示了每個(gè)基因模塊與OA的相關(guān)顯著性，其中棕色基因模塊與OA的相關(guān)關(guān)系最為顯著。圖3C顯示了棕色模塊基因與OA相關(guān)的顯著性，因此將棕色基因模塊作為關(guān)鍵基因模塊。

2.3 關(guān)鍵基因模塊的功能富集分析GO富集分析結(jié)果表明，棕色基因模塊的基因主要富集在脂類代謝調(diào)控、T細(xì)胞活化、細(xì)胞凋亡、炎癥反應(yīng)調(diào)控、組織缺氧反應(yīng)、細(xì)胞饑餓反應(yīng)等生物學(xué)過(guò)程(圖4)。KEGG通路富集分析結(jié)果表明，棕色模塊基因主要富集在MAPK信號(hào)通路、破骨細(xì)胞通路、腫瘤信號(hào)通路、脂肪細(xì)胞因子信號(hào)通路、mTOR信號(hào)通路、晝夜節(jié)律等信號(hào)通路(圖5)。

2.4 關(guān)鍵基因模塊的PPI網(wǎng)絡(luò)分析采用STRING在線數(shù)據(jù)庫(kù)以及Cytoscape軟件構(gòu)建棕色模塊基因的PPI(圖6A)，利用cytoHubba插件的degree算法篩選出PPI前10個(gè)基因(圖6B)，分別是JUN原癌基因(jun proto-oncogene, JUN)、CCAAT增強(qiáng)子結(jié)合蛋白β(CCAAT enhancer binding protein beta，CEBPB)、Fos原癌基因(Fos proto-oncogene，F(xiàn)OS)、RELA原癌基因(RELA proto-oncogene，RELA)、泛素C(Ubiquitin C，UBC)、CCAAT增強(qiáng)子結(jié)合蛋白α(CCAAT enhancer binding protein alpha，CEBPA)、激活轉(zhuǎn)錄因子3(activating transcription factor 3，ATF3)、DNA 損傷誘導(dǎo)的轉(zhuǎn)錄因子3(DNA damage inducible transcript 3，DDIT3)、BCL6轉(zhuǎn)錄抑制因子(BCL6 transcription repressor，BCL6)和干擾素調(diào)節(jié)因子4(interferon regulatory factor 4，IRF4)。

2.5 關(guān)鍵基因的表達(dá)水平在cytoHubba插件的degree算法篩選出PPI前10個(gè)基因中，DDIT3和BCL6基因目前未見(jiàn)其與骨關(guān)節(jié)炎相關(guān)的報(bào)道研究，因此最終選取DDIT3和BCL6作為關(guān)鍵基因。圖7顯示了DDIT3和BCL6在骨關(guān)節(jié)炎與正常對(duì)照兩組間的表達(dá)量水平情況，這兩個(gè)關(guān)鍵基因在骨關(guān)節(jié)炎組樣本中相對(duì)于正常對(duì)照組表達(dá)量顯著下調(diào)。

2.6 關(guān)鍵基因的驗(yàn)證在驗(yàn)證集GSE57218中，DDIT3和BCL6這兩個(gè)關(guān)鍵基因的表達(dá)量如圖8A、B所示，DDIT3和BCL6的表達(dá)量在OA組中相對(duì)于正常對(duì)照組顯著下調(diào)。在另一個(gè)驗(yàn)證集GSE169077中，DDIT3和BCL6的表達(dá)量如圖8C、D所示，其表達(dá)量在OA組中相對(duì)于正常對(duì)照組也顯著下調(diào)。

3 討 論

骨關(guān)節(jié)炎是一種退行性骨關(guān)節(jié)疾病，以軟骨退化和破壞、關(guān)節(jié)滑膜炎、軟骨下骨改變及骨贅形成為病理特征[1]，導(dǎo)致日常生活能力嚴(yán)重下降，給患者和醫(yī)療保健系統(tǒng)帶來(lái)沉重的經(jīng)濟(jì)負(fù)擔(dān)[2]。目前，骨關(guān)節(jié)炎的發(fā)病機(jī)制仍不清楚，了解其發(fā)病機(jī)制對(duì)骨關(guān)節(jié)炎的診斷、治療和預(yù)后具有重要的意義[18-19]，本研究從GEO數(shù)據(jù)庫(kù)下載GSE114007的轉(zhuǎn)錄組數(shù)據(jù)，利用“DESeq2”包進(jìn)行差異分析，篩選出1752個(gè)DEGs，其中上調(diào)基因927個(gè)，下調(diào)基因825個(gè)。WGCNA分析方法將DEGs構(gòu)建基因共表達(dá)網(wǎng)絡(luò)，得到6個(gè)不同表達(dá)模式的基因模塊，其中棕色基因模塊(共281個(gè)基因)與OA的相關(guān)系數(shù)最高，且呈負(fù)相關(guān)，表明了棕色基因模塊是參與骨關(guān)節(jié)炎發(fā)生發(fā)展的關(guān)鍵基因模塊。GO富集分析表明，棕色模塊基因主要富集在脂類代謝調(diào)控、T細(xì)胞活化、細(xì)胞凋亡、炎癥反應(yīng)和組織缺氧反應(yīng)等生物調(diào)控過(guò)程，在這些生物學(xué)過(guò)程中，脂類代謝、細(xì)胞凋亡和炎癥反應(yīng)已經(jīng)被證明與OA的發(fā)展有關(guān)[20-23]。KEGG通路富集分析表明，MAPK信號(hào)通路、破骨細(xì)胞通路、腫瘤信號(hào)通路、脂肪細(xì)胞因子信號(hào)通路和mTOR信號(hào)通路與OA的發(fā)病機(jī)制有關(guān)，而這些信號(hào)通路已經(jīng)被證明在OA發(fā)病機(jī)制中起著重要的作用[24-25]。正常軟骨組織中的代謝受嚴(yán)格的調(diào)節(jié)，以維持合成和降解之間的平衡，有研究認(rèn)為骨關(guān)節(jié)炎發(fā)生時(shí)，這種平衡被打破，在骨關(guān)節(jié)炎的軟骨退化后其代謝和基因表達(dá)譜都會(huì)發(fā)生變化。

PPI網(wǎng)絡(luò)分析篩選出互作關(guān)系最強(qiáng)的10個(gè)基因JUN、CEBPB、FOS、RELA、UBC、CEBPA、ATF3、DDIT3、BCL6、IRF4。有研究發(fā)現(xiàn)JUN通過(guò)調(diào)節(jié)軟骨細(xì)胞合成代謝和分解代謝基因的表達(dá)在OA軟骨破壞中發(fā)揮重要作用[26]。LU等人[27]進(jìn)一步研究發(fā)現(xiàn)JUN參與軟骨細(xì)胞的凋亡，且與OA的發(fā)生有緊密關(guān)系。FOS是一種具有調(diào)控細(xì)胞的增殖、分化和轉(zhuǎn)化生理功能的原癌基因[28]，其可上調(diào)MMP1、抑制TIMP1啟動(dòng)子的轉(zhuǎn)錄，從而降解軟骨組織，是其致骨關(guān)節(jié)炎發(fā)生的機(jī)制之一[29]。ATF3與多種疾病的發(fā)病機(jī)制有關(guān)，包括代謝性、免疫性和炎癥性疾病，研究發(fā)現(xiàn)ATF3通過(guò)核因子-kB(NF-kB)途徑軟骨細(xì)胞炎性細(xì)胞因子的表達(dá)參與OA的發(fā)病過(guò)程[30]。此外，CEBPB、RELA、UBC、CEBPA和IRF4這些基因也陸續(xù)被發(fā)現(xiàn)與OA的發(fā)病機(jī)制有關(guān)[31-33]。

DDIT3是一種內(nèi)質(zhì)網(wǎng)應(yīng)激通路中的關(guān)鍵轉(zhuǎn)錄因子，研究發(fā)現(xiàn)DDIT3通過(guò)內(nèi)質(zhì)網(wǎng)應(yīng)激通路調(diào)控細(xì)胞凋亡和炎癥反應(yīng)過(guò)程[34]，XIONG等人[35]進(jìn)一步研究發(fā)現(xiàn)在軟骨發(fā)育過(guò)程中DDIT3通過(guò)內(nèi)質(zhì)網(wǎng)應(yīng)激通路參與細(xì)胞凋亡，以維持正常的代謝。YU等人[36]研究發(fā)現(xiàn)DDIT3通過(guò)誘導(dǎo)SOX9的表達(dá)來(lái)調(diào)控軟骨細(xì)胞分化，YANG等人[37]進(jìn)一步研究闡明了DDIT3通過(guò)SIRT1-AKT途徑促進(jìn)軟骨細(xì)胞的自噬，這些研究表明DDIT3對(duì)軟骨的代謝、分化和凋亡具有極其重要的作用。GO和KEGG富集分析顯示DDIT3參與脂類代謝、細(xì)胞凋亡和MAPK信號(hào)通路，因此DDIT3有可能通過(guò)脂類代謝、細(xì)胞凋亡和MAPK信號(hào)通路參與OA的發(fā)生發(fā)展。此外，PPI網(wǎng)絡(luò)分析表明，DDTI3與JUN、FOS、ATF3存在著相互作用網(wǎng)絡(luò)關(guān)系，因此DDIT3與這些基因存在相互調(diào)節(jié)的作用，又因JUN、FOS、ATF3基因與OA的發(fā)病機(jī)制有關(guān)，所以DDIT3可能通過(guò)這些基因進(jìn)而參與OA的發(fā)病機(jī)制。有研究發(fā)現(xiàn)FOS和DDIT3存在相互調(diào)控關(guān)系共參與軟骨形成和分化[38]，LI等[39]研究發(fā)現(xiàn)DDIT3和JUN通過(guò)形成新的復(fù)合物參與細(xì)胞凋亡的調(diào)控。此外，DDIT3通過(guò)IL-1β和MMP3介導(dǎo)軟骨細(xì)胞分解代謝和細(xì)胞凋亡反應(yīng)[40]。

BCL6是一種轉(zhuǎn)錄抑制因子，在細(xì)胞增殖、炎癥反應(yīng)和凋亡中起到重要的作用。GO富集分析顯示BCL6參與炎癥反應(yīng)、免疫反應(yīng)、衰老和細(xì)胞凋亡等生物學(xué)過(guò)程。有研究發(fā)現(xiàn)BCL6對(duì)成骨細(xì)胞分化有重要的作用[41]，另外，有研究發(fā)現(xiàn)BCL6對(duì)破骨細(xì)胞分化也起到重要的作用[42]。許多研究發(fā)現(xiàn)BCL6在類風(fēng)濕性關(guān)節(jié)炎發(fā)病機(jī)制中扮演著重要的角色[43]。PPI網(wǎng)絡(luò)分析顯示BCL6與JUN、FOS及IRF4存在蛋白互作關(guān)系，因此BCL6有可能通過(guò)JUN、FOS和IRF4之間相互作用參與OA的發(fā)病機(jī)制。CHAWEEWANNAKORN等[44]研究發(fā)現(xiàn)BCL6通過(guò)下調(diào)FOS影響破骨細(xì)胞調(diào)節(jié)因子的生成。這些研究表明BCL6與骨關(guān)節(jié)炎的發(fā)生有著密切關(guān)聯(lián)。

本研究又在GSE57218和GSE169077兩個(gè)數(shù)據(jù)集中進(jìn)行驗(yàn)證，DDIT3和BCL6在兩個(gè)驗(yàn)證集的骨關(guān)節(jié)炎組中的表達(dá)量也是顯著低于正常對(duì)照組，進(jìn)一步驗(yàn)證了本研究的結(jié)果，因此，我們可以到得DDIT3和BCL6與骨關(guān)節(jié)炎的發(fā)生發(fā)展有密切關(guān)聯(lián)。雖然本研究結(jié)果在另外兩個(gè)表達(dá)譜數(shù)據(jù)得到了驗(yàn)證，但我們研究還有不足之處。其一DDIT3和BCL6這兩個(gè)關(guān)鍵基因沒(méi)有進(jìn)一步實(shí)驗(yàn)驗(yàn)證，其二是對(duì)DDIT3和BCL6關(guān)鍵基因參與OA發(fā)病的具體機(jī)制沒(méi)有進(jìn)一步地研究分析，只是通過(guò)GO、KEGG富集分析和現(xiàn)有文獻(xiàn)研究結(jié)果進(jìn)行關(guān)聯(lián)探尋其參與OA的潛在發(fā)病機(jī)制。

綜上所述，本研究結(jié)果表明DDIT3和BCL6與骨關(guān)節(jié)炎的發(fā)生發(fā)展有密切關(guān)聯(lián)，DDIT3和BCL6可能通過(guò)MAPK信號(hào)通路、mTOR信號(hào)通路、脂類代謝參與骨關(guān)節(jié)炎的發(fā)生發(fā)展，因此DDIT3和BCL6可作為新的骨關(guān)節(jié)炎候選生物標(biāo)志物。