基于WGCNA篩選骨關節(jié)炎的生物標志物

龍通華，張紅參，許冬梅，曾志華，傅子原，關皓鄴

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

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

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

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

1 材料與方法

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

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

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

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

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

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

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

2 結 果

2.1 差異表達基因篩選在骨關節(jié)炎和正常對照兩組中共篩選出1752個DEGs，其中上調基因927個，下調基因825個，DEGs的結果展示如圖1。

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

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

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

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

2.6 關鍵基因的驗證在驗證集GSE57218中，DDIT3和BCL6這兩個關鍵基因的表達量如圖8A、B所示，DDIT3和BCL6的表達量在OA組中相對于正常對照組顯著下調。在另一個驗證集GSE169077中，DDIT3和BCL6的表達量如圖8C、D所示，其表達量在OA組中相對于正常對照組也顯著下調。

3 討 論

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

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

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

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

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

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