基于WGCNA篩選骨關節炎的生物標志物

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

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

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

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

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

1 材料與方法

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

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

1.3 構建DEGs共表達基因模塊利用“WGCNA”包[5]對DEGs構建共表達網絡，首先pickSoftThreshold函數計算出相關系數為0.86時對應的軟閾值為10。以軟閾值為10將差異表達基因構建成鄰接矩陣，之后轉變為拓撲重疊矩陣[14]，以識別基因模塊。基因模塊的最小基因數為30，合并相似基因模塊的切割高度的閾值設為0.25[5]。

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

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

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

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

2 結 果

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

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

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

2.4 關鍵基因模塊的PPI網絡分析采用STRING在線數據庫以及Cytoscape軟件構建棕色模塊基因的PPI(圖6A)，利用cytoHubba插件的degree算法篩選出PPI前10個基因(圖6B)，分別是JUN原癌基因(jun proto-oncogene, JUN)、CCAAT增強子結合蛋白β(CCAAT enhancer binding protein beta，CEBPB)、Fos原癌基因(Fos proto-oncogene，FOS)、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)和干擾素調節因子4(interferon regulatory factor 4，IRF4)。

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

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

3 討 論

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

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

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

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

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

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