生物信息學分析篩選骨關節炎關鍵基因的研究▲

蔣 捷 黃林科，4 韋林華，5 徐家科，3* 陳蔚蔚#

(1 廣西醫科大學第一附屬醫院創傷骨科手外科，廣西南寧市 530021；2 廣西醫科大學再生醫學研究中心，廣西南寧市 530021；3 西澳大利亞大學生物醫學與科學學院，澳大利亞珀斯 6009；4 廣西醫科大學第二附屬醫院骨科，廣西南寧市 530006；5 廣西醫科大學附屬南寧市傳染病醫院外科，廣西南寧市 530023)

骨關節炎(osteoarthritis，OA)是影響中老年人生活質量最常見的慢性病之一，發病率為6%～32%[1]。OA是一種異質性疾病，具體病因不詳，一般認為與遺傳、年齡、肥胖及機械應力等有關[2]。OA病理變化主要包括4個方面：關節軟骨退變、軟骨下骨硬化、關節周圍骨贅形成和滑膜炎癥浸潤[3]。對OA發病機制和治療手段的研究，目前主要集中在軟骨組織和軟骨細胞上[4]，但研究結果并不理想。近年來，通過對滑膜組織及軟骨下骨在OA發生發展過程中的作用的研究，得到了更多的認識[5]。滑膜炎引起的不正常滑膜液循環導致關節積液、引起關節疼痛，是早期OA的主要表現[6]。不同階段OA伴有不同程度的滑膜炎癥[7]。因此，滑膜炎癥的研究對揭示OA的發病機制、疾病進展和治療方案的選擇有著重要的作用。

隨著生物信息學的快速發展，越來越多的研究表明OA的發生發展與一系列的基因和信號通路相關[8]。利用GEO數據庫(Gene Express Omnibus, http: //www.ncbi.nlm.nih.gov/geo)中關節滑膜基因表達譜芯片數據，篩選健康群體和OA患者滑膜組織的差異表達基因(differentially expressed genes，DEGs)，對DEGs進行聚類分析及功能富集分析，并進一步構建蛋白相互作用(protein-protein interaction, PPI)網絡，篩選出與OA密切相關的關鍵基因，旨在進一步識別、鑒定OA發生發展中的靶基因，為OA發病分子機制的明確及分子標志物的篩選提供理論指導。現報告如下。

1 材料與方法

1.1 基因芯片數據 從GEO數據庫下載基于GPL96[HG-U133A] Affymetrix Human Genome U133A Array平臺的基因表達譜芯片數據集：GSE55457、GSE55235、GSE12021。3個數據集共包括了30例OA患者和29例健康對照者的滑膜組織(見表1)。所有數據集中的類風濕關節炎患者樣本被排除，GSE12021數據集中基于GPL97平臺的樣本也被排除。

表1 3個GEO數據集信息

1.2 DEGs分析 利用GEO自帶的在線工具GEO2R(https://www.ncbi.nlm. nih.gov/geo/geo2r/)對每個數據集的OA患者及健康對照者的滑膜組織進行DEGs分析，按照調整后P值(adjustedP)<0. 05，|log2FC|>1(FC表示DEGs上調或下調倍數)的篩選條件分別篩選3個數據集中OA和健康對照者滑膜組織的DEGs。使用在線分析工具jvenn(http://jvenn.toulouse.inra.fr/app/index.html)對3個數據集中共同的DEGs進行篩選，確定DEGs。

1.3 DEGs的GO功能富集及KEGG通路富集分析 利用在線分析工具DAVID 6.8(https://david.ncifcrf.gov/)[9]對篩選的DEGs進行GO(Gene Ontology, GO)功能富集分析和KEGG(Kyoto Encyclopedia of Genes and Genomes,KEGG)通路富集分析，以P<0.05為差異有統計學意義。

1.4 差異基因PPI網絡的構建與關鍵基因分析 應用在線分析數據庫STRING(http://string-db.org/cgi/input.pl)[10]，將篩選出的DEGs導入數據庫，進行差異基因PPI分析。利用Cytoscape軟件構建PPI網絡[11]，并應用其中的Cytohubba插件選取前10 位的差異基因作為蛋白質網絡中的hub基因。MCODE插件對PPI進行模塊化分析(參數：degree cutoff=2, node score cutoff=0.2, K-core=2, max. cepth=100)，選取score>3且基因數大于4個的子網絡作為顯著性模塊。

2 結 果

2.1 健康對照者與OA患者滑膜組織的DEGs 分別從基因芯片數據集GSE55457、GSE55235、GSE12021中篩選出786、1 048、808個差異基因。在3個數據集中均有的DEGs為126個(圖1)。

圖1 GSE55457、GSE55235、GSE12021的DEGs

2.2 DEGs的GO功能富集和KEGG通路富集分析結果 GO富集分析表明，滑膜組織中DEGs參與的主要生物學過程有：細胞對促腎上腺皮質激素釋放激素刺激的反應、通過mRNA剪接體對mRNA的剪接、負向調控基因表達、正向調控單核細胞聚集、正向調控成纖維細胞增殖、蛋白轉運、調控炎癥反應、骨骼肌細胞分化、調節呼吸過程中氣體交換、正向調控細胞遷移；細胞成分包括核質、細胞核、生長錐膜；分子功能主要集中于核苷酸結合、肝素結合、RNA聚合酶Ⅱ轉錄因子活性與配體激活序列特異性DNA結合、類固醇激素受體活性。詳見表2。

表2 DEGs的GO分析結果

KEGG富集通路分析顯示，滑膜組織DEGs的主要信號通路富集在包括MAPK信號通路、EB(Epstein-Barr，EB)病毒分子感染通路、Ⅰ型人類T細胞白血病病毒分子感染通路、癌癥中蛋白聚糖分子通路、惡性腫瘤分子通路、破骨細胞分化通路、PI3K/Akt信號通路、胰島素信號轉導通路、腫瘤中miRNAs分子通路、ErbB信號通路、HIF-1信號通路和膀胱癌分子相關通路。詳見圖2。

圖2 DEGs參與的KEGG信號通路

2.3 DEGs的PPI網絡分析 將DEGs導入STRING和Cytoscape構建出PPI網絡，該網絡由84個結點和233條邊構成(圖3)。Cytohubba插件確認的度值前10的hub基因是(圖4)：JUN、BTG2、ATF3、DUSP1、HNRNPA1、SF1、SRSF7、TRA2B、TRA2A、PPP1R15A。

圖3 DEGs編碼PPI網絡

圖4 hub基因及其擴展PPI網絡

MCODE插件分析出score>3且基因數多于4個的子網絡有4個(圖5)，其中最顯著的相互作用模塊(score:7)由7個結點和21條邊構成，包含了TRA2B、SRRM2、SRSF11、TRA2A、HNRNPA1、SRSF7、SF1 7個hub基因，這些基因存在非常復雜的相互作用。

圖5 MCODE確認的網絡核心模塊

3 討 論

OA是世界范圍內致殘的主要因素[12]。但OA發病的病理機制仍然不明確，因此了解OA的發病機制將有助于早期診斷、判斷預后和發現藥物的新靶點[13]。

本研究從GEO數據庫獲取健康對照者和OA患者滑膜組織的基因表達譜芯片數據，利用一系列在線分析工具(GEO2R、jvenn、DAVID、STRING、Cytoscape)，結合生物信息學方法對二者的DEGs進行分析，篩選出了126個DEGs。對DEGs進行GO功能富集，這些基因主要參與細胞對促腎上腺皮質激素釋放激素刺激的反應、通過mRNA剪接體對mRNA的剪接、負向調控基因表達、正向調控單核細胞聚集、正向調控成纖維細胞增殖等生物學過程。通過KEGG通路富集分析，我們確定了與本團隊前期研究相一致的影響OA發病機制通路，包括MAPK、PI3K/Akt和HIF-1信號通路。研究表明，抑制MAPK信號通路可以減輕滑膜炎癥，改善炎癥誘導的軟骨細胞凋亡，促進軟骨再生[14-15]。此外，PI3K/Akt信號通路對細胞的增殖、凋亡過程具有重要的調節作用，在OA的發生、發展過程中起著一定的作用。有研究報道PI3K/Akt信號通路在OA中高表達[16]，抑制PI3K/Akt信號通路會降低膠原誘導關節炎小鼠滑膜鈣粘蛋白-11的表達和大鼠關節炎的炎癥反應[17-18]。HIF-1是細胞對缺氧反應的主要轉錄調節因子，由于軟骨組織生長在一個低氧的微環境中，HIF-1在軟骨穩態的平衡維持中發揮重要作用，活化的HIF-1與軟骨下骨重構、新生血管生成有密切關系[19-20]。

同樣，我們構建了DEGs的PPI網絡，利用Cytoscape軟件中的Cytohubba插件對差異基因PPI網絡中連接最為緊密且位于重要結點位置的基因進行分組，篩選出了JUN、BTG2、ATF3、DUSP1、HNRNPA1、SF1、SRSF7、TRA2B、TRA2A、PPP1R15A等hub基因。JUN是編碼激活蛋白-1的亞單位，白細胞介素-1能夠通過激活JUN誘導基質金屬蛋白酶-13產生，從而導致軟骨基質降解[21-22]。JUN編碼的c-JUN蛋白對調節軟骨細胞凋亡起著重要的作用[23]。ATF3在OA軟骨中的表達顯著升高，ATF3可以調節炎癥因子如白介素細胞-6的分泌，影響軟骨的代謝[24]。TRA2B是一種重要的RNA結合蛋白，在預后較差的骨肉瘤患者中高表達，miR-206可下調TRA2B表達誘導骨肉瘤細胞凋亡[25]。BTG2屬于BTG/TOB家族成員，可有效下調PI3K和Akt磷酸化抑制骨肉瘤細胞遷移侵襲能力[26]。敲除BTG2導致小鼠的腰椎后移，與BMP信號通路活化減弱相關[27-28]。HNRNPA1是一種RNA結合蛋白，參與信使RNA的代謝和轉運，突變可能引起多系統蛋白病。TRAF6調節HNRNPA1的泛素化失衡會誘發骨髓異常增生綜合征[29-30]。SRSF7、SRRM2、SRSF11、BTG2、SF1、HNRNPA1、PPP1R15A暫時未見與OA相關的報道，可能是新的潛在的OA發生相關基因。

在本研究中，我們同時發現了一些與OA有關的新的標志物。DUSP1在3個數據集的OA患者中相對健康對照者表達量明顯減少。DUSP1是絲裂原活化蛋白激酶(mitogen-activated protein kinases, MAPKs)的內源性抑制物，下調或者敲除DUSP1可以增強MAPKs作用[31-32]。在OA中，白細胞介素-1β誘導的MAPKs的活化導致炎癥因子誘導的軟骨細胞凋亡，基質金屬蛋白酶-3、基質金屬蛋白酶-13等酶的大量產生導致軟骨基質降解[33-34]。因此，是否可以通過增加DUSP1的表達，減少MAPKs的活化，從而達到減輕炎癥、保護軟骨、延緩OA發生的作用，值得進一步研究。

本研究的不足之處：(1)未能進一步做到hub基因的GO和KEGG分析；(2)未在實例標本中證明hub基因的表達情況；(3)未能分析MCODE插件中得出的各子網絡的作用。在今后的研究中，我們會完善相關實驗，得出更完整的結論。

綜上所述，在對本研究篩選得到的DEGs如JUN、ATF3、DUSP1等這些分子生物功能及信號通路的進一步研究中，對OA的診治將會有更多的新發現。BTG2、HNRNPA1、SF1、SRSF7、TRA2B、TRA2A、PPP1R15A這些基因可能是潛在的影響OA的發生發展的新基因，雖然其在OA中的價值仍然需要通過大量的實驗證實，但為今后的進一步研究提供了重要的線索和參考信息。