基于生物信息學的骨關(guān)節(jié)炎差異基因譜預(yù)測

祝乃強 侯靜怡 劉金欣 馬桂云

(1承德醫(yī)學院附屬醫(yī)院脊柱外二科，河北 承德 067000;2承德醫(yī)學院)

骨關(guān)節(jié)炎(OA)多見于中老年人，65歲以上老年人患病率高達50%〔1〕。近年來，研究發(fā)現(xiàn)OA發(fā)病機制并不是一組靶點或一些基因間的相互作用關(guān)系出現(xiàn)了明顯改變，而是一組靶點或一些基因間的相互作用關(guān)系出現(xiàn)了明顯改變從而導(dǎo)致其相關(guān)功能模塊或生物網(wǎng)絡(luò)發(fā)生故障，其中包括軟骨細胞外基質(zhì)的酶降解，新細胞外基質(zhì)合成缺乏，細胞死亡與凋亡，軟骨細胞的異?；罨c分化等〔2，3〕。因此，尋找新的更有意義的生物標志物及其參加的信號通路，在預(yù)防與治療OA疾病發(fā)生和發(fā)展過程中起到重要作用。隨著高通量技術(shù)的快速發(fā)展，生物數(shù)據(jù)呈現(xiàn)指數(shù)級的增長，以基因芯片為代表的高通量生物檢測技術(shù)得到廣泛應(yīng)用，其可在短時間內(nèi)獲得大量的疾病相關(guān)基因信息，使得疾病在基因?qū)用孢M行全面而細致的分析成為可能〔4〕。本研究采用生物信息學方法分析OA基因芯片數(shù)據(jù)，目的在于尋找OA病理過程中的關(guān)鍵基因及其參與信號通路，為OA早期診斷、早期治療干預(yù)、藥物研發(fā)提供新的思路。

1 材料與方法

1.1數(shù)據(jù)收集與預(yù)處理 以“Osteoarthritis”作為檢索詞，在NCBI的基因表達匯編(GEO)數(shù)據(jù)庫檢索相應(yīng)的基因芯片，下載符合實驗要求的4套數(shù)據(jù)集，其系列號為GSE12021、GSE55235、GSE55457和GSE55584，均來自GPL96平臺，芯片類型為Affymetrix Human Genome U133a芯片。GSE12021芯片數(shù)據(jù)集由Huber等〔5〕于2008年提交，該數(shù)據(jù)集包括13例正?；そM織樣本，20例OA組織樣本，20例類風濕性關(guān)節(jié)炎滑膜組織樣本；GSE55235、GSE55457和GSE55584數(shù)據(jù)集是由Woetzel等〔6〕于2014年提交的，共包含20例正?；そM織樣本，26例OA滑膜組織樣本，33例類風濕性關(guān)節(jié)炎滑膜組織樣本。本研究共選取4個數(shù)據(jù)集中的正?；そM織樣本和OA滑膜組織樣本作為后續(xù)分析的樣本。

1.2差異表達基因的篩選分析 利用R語言中的 “sva”分析包〔7〕對GSE12021、GSE55235、GSE55457和GSE55584數(shù)據(jù)集中的數(shù)據(jù)進行合并及批次校正，利用“l(fā)imma”分析包〔8〕對上述4個數(shù)據(jù)集進行分析，以P<0.05，|log2FC|>1為篩選條件選擇數(shù)據(jù)集的差異表達基因。

1.3差異基因表達的基因本體論(GO)功能富集和京都基因與基因組百科全書(KEGG)通路富集分析 GO分析用于高通量基因組或轉(zhuǎn)錄組數(shù)據(jù)鑒定的生物學機制，包括分子功能(MF)、生物學過程(BP)、細胞組分(CC)〔9〕；KEGG 數(shù)據(jù)庫用于識別候選目標的功能及生物相關(guān)性〔10〕。本研究利用R語言中的“clusterProfiler”富集分析包〔11〕對差異基因進行GO功能富集分析和KEGG通路富集分析。

1.4蛋白-蛋白互作網(wǎng)絡(luò)(PPI)構(gòu)建及關(guān)鍵基因分析(PPI) 將所有的差異基因?qū)氲絊TRING線上數(shù)據(jù)庫(https：//string-db.org/)〔12〕，進行蛋白相互作用分析，并利用Cytoscape軟件對PPI進行可視化。應(yīng)用Cytoscape軟件中的Cytohubba插件對網(wǎng)絡(luò)進行拓撲學分析，以度值(degree)前10位的差異基因作為PPI中的關(guān)鍵基因。DisGeNET數(shù)據(jù)庫(http：//www.disgenet.org/web/DisGeNET/me，Version 5.0)〔13〕是包含與人類疾病相關(guān)的基因和變體的平臺之一。在 DisGeNET 數(shù)據(jù)庫中選用基因進行檢索，將上述靶點依次導(dǎo)入到數(shù)據(jù)庫中，獲取靶點類型信息(Protein Class)。

1.5關(guān)鍵模塊分析及功能富集分析 利用Cytoscape軟件中的插件復(fù)雜分子檢測(MCODE)對PPI進行模塊分析，標準設(shè)置：節(jié)點得分截止值=0.2，K核心值(K-Core)=2，度值截止=2，計算MCODE的得分值；以MCODE分數(shù)>6分作為顯著性模塊的篩選標準，并且利用DAVID(https：//david.ncifcrf.gov/)在線數(shù)據(jù)庫對篩選出的模塊中的基因進行通路富集分析〔14〕。

2 結(jié) 果

2.1差異表達基因的篩選結(jié)果 GSE12021、GSE55235、GSE55457和GSE55584的數(shù)據(jù)集共包含33例正常滑膜組織樣本和40例OA滑膜組織樣本，篩選出符合條件的差異基因共320個，其中上調(diào)基因172個，下調(diào)基因148個。

2.2差異基因的GO功能富集和KEGG通路富集分析 將差異基因進行GO富集分析，其中上調(diào)的差異基因共有151個條目富集，主要包括加雙氧酶活動、鐵離子結(jié)合、細胞黏合力、泛素蛋白連接酶集合等，；下調(diào)的差異基因共有216個條目富集，主要包括抑制轉(zhuǎn)運因子結(jié)合、核苷酸糖跨膜轉(zhuǎn)運活動，激素受體結(jié)合等，詳細內(nèi)容見圖1、圖2。

圖1 上調(diào)差異基因的GO功能富集分析

圖2 下調(diào)差異基因的GO功能富集分析

KEGG通路富集結(jié)果表明，上調(diào)的差異基因共有159個條目富集，主要包括絲裂原活化蛋白激酶(MAPK)信號通路、脂肪酸代謝通路、凋亡通路等；下調(diào)的基因共有133個條目富集，主要包括Janus激酶/信號轉(zhuǎn)導(dǎo)與轉(zhuǎn)錄激活子(JAK-STAT)信號通路、原發(fā)性免疫通路、造血細胞譜系通路等，詳細內(nèi)容見圖3、4。

圖3 上調(diào)差異基因的KEGG通路富集分析

2.3PPI及關(guān)鍵基因分析 將篩選出的差異基因?qū)氲絊TRING在線數(shù)據(jù)庫中，得到差異基因的PPI，該網(wǎng)絡(luò)包含313個節(jié)點，2 160條邊。如圖5所示，該網(wǎng)絡(luò)包含167個上調(diào)基因和144個下調(diào)基因。應(yīng)用Cytoscape軟件的Cytohubba插件，分析網(wǎng)絡(luò)的拓撲學參數(shù)，且以degree值的大小表示該基因在此網(wǎng)絡(luò)中的重要性，如若degree越大，表明與其關(guān)聯(lián)的基因越多，則其基因越重要。

圖4 下調(diào)差異基因的KEGG通路富集分析

(白色節(jié)點代表上調(diào)基因，灰色節(jié)點代表下調(diào)基因)圖5 差異基因的蛋白-蛋白互相作用網(wǎng)絡(luò)

以degree前10位(degree > 40)的基因作為網(wǎng)絡(luò)的關(guān)鍵基因，所篩選的關(guān)鍵基因為：連環(huán)蛋白(CTNNB)1，糖原合成酶激酶(GSK)-3β，信號轉(zhuǎn)導(dǎo)和轉(zhuǎn)錄激活因子(STAT)1，Ras同系物集組成員(RHO)C，組蛋白去乙酰化酶(HDAC)9，早老素(PSEN)1，賴氨酸脫甲基酶(KDM)5C，β-分泌酶(BACE)1，JAK 3和卡林蛋白(CUL1)；所涉及的OA過程中有信號通路、轉(zhuǎn)錄因子、蛋白(鈣結(jié)合蛋白)、酶(蛋白酶、激酶、酶調(diào)節(jié)劑、連接酶)。

2.4功能模塊分析 應(yīng)用Cytoscape插件MCODE尋找差異基因PPI的關(guān)鍵模塊，其中符合MCODE得分大于4分的模塊共有3個，模塊1的MCODE得分為7.867，16個節(jié)點和59條邊。見圖6。模塊1中的基因進行功能富集分析，模塊1中基因主要參與6條富集通路，其中BP 3個條目，包括轉(zhuǎn)錄/DNA模板化、神經(jīng)元分化負調(diào)劑、RNA聚合酶Ⅱ啟動子的轉(zhuǎn)錄正調(diào)節(jié)；CC 1個條目-質(zhì)膜外側(cè)和MF 2個條目，包括RNA聚合酶Ⅱ核心啟動子序列、轉(zhuǎn)錄因子活性-DNA結(jié)合；KEGG通路富集分析顯示，模塊1中的基因主要參與5條信號通路，如JAK-STAT信號通路、趨化因子、磷脂酰肌醇3激酶(PI3K)-絲氨酸-蘇氨酸激酶(Akt)信號通路等。見表1。

(淺色節(jié)點代表上調(diào)基因，深色節(jié)點代表下調(diào)基因)圖6 功能模塊的PPI

表1 功能模塊的功能富集和通路富集分析

FDR：錯誤發(fā)現(xiàn)率

3 討 論

OA是最常見的一種全關(guān)節(jié)疾病，肥胖、年齡、過度運動、炎癥、遺傳及創(chuàng)傷等多種因素與OA的發(fā)生密切聯(lián)系〔15〕。目前，OA的病因及發(fā)病機制尚不明確。分析OA的潛在發(fā)病機制對于診斷預(yù)后和鑒定藥物治療的靶點至關(guān)重要。由于高通量測序和微陣列技術(shù)可同時提供有關(guān)人類基因組中數(shù)千種基因表達水平的信息，因此該方法已廣泛用于預(yù)測OA的潛在診斷和治療靶點。

本研究進行差異基因分析表明OA發(fā)病機制的復(fù)雜性，其并不是簡單的單個基因/靶點單獨其作用，而是多基因/靶點相互作用的復(fù)雜結(jié)果。KEGG通路富集結(jié)果顯示，這些差異基因主要涉及MAPK信號通路，脂肪酸代謝通路，凋亡通路、JAK-STAT信號通路，原發(fā)性免疫通路，造血細胞譜系通路等。MAPK信號通路能夠調(diào)控軟骨細胞的增殖、凋亡、胞外基質(zhì)代謝、炎癥因子分泌等過程，在OA病理過程中發(fā)揮重要作用〔16，17〕。Ge等〔16〕研究表明，3-磷酸肌醇依賴性蛋白激酶(PDK)1可通過p38 MAPK信號通路促進OA的軟骨細胞凋亡。而Peng等〔18〕研究進一步發(fā)現(xiàn)OA中的雙特異性磷酸酶(DUSP)1過表達抑制了MAPK信號通路的激活和OA相關(guān)基質(zhì)的表達。研究表明，JAK/STAT信號通路受到促炎細胞因子基因表達而顯著上調(diào)的調(diào)控。Lim等〔19〕研究表明，MMP-13在IL-1β處理的人軟骨肉瘤細胞系中的表達誘導(dǎo)了JAK2和STAT1/STAT2激活及MMP-13基因表達，其被泛酪氨酸激酶抑制劑AG490阻斷。然而，也有研究發(fā)現(xiàn)〔20〕，用IL-1β處理這種軟骨細胞系也激活了p38-MAPK/c-Fos/AP-1。

利用STRING在線數(shù)據(jù)庫對這些差異基因進行PPI網(wǎng)絡(luò)分析得到位于中心的關(guān)鍵靶點，其中，STAT1作為調(diào)節(jié)巨噬細胞炎癥反應(yīng)中不可或缺的組成部分，如若激活的STAT1被磷酸并易位到細胞核中以上調(diào)炎癥的基因〔21〕。JAK3作為一種非受體酪氨酸激酶，可與細胞因子受體復(fù)合物的常見γ鏈(IL2RG)相互作用〔22〕。

PI3K-Akt信號通路對于生理和病理條件下細胞生長和存活的許多方面都至關(guān)重要〔23〕。研究表明，細胞內(nèi)信號轉(zhuǎn)導(dǎo)途徑MAPK和PI3K-Akt已被證實在調(diào)節(jié)炎癥介質(zhì)和骨關(guān)節(jié)炎發(fā)病機制中的細胞外基質(zhì)改變中都是至關(guān)重要的〔24，25〕。其中IL-1β能夠激活MAPK的三種途徑(ERK，p38，JNK)及PI3K/Akt信號通路，隨著誘導(dǎo)分解代謝的表達，可進一步促進細胞外基質(zhì)破壞〔26〕，表明MAPK和PI3K/Akt的調(diào)節(jié)在IL-1β處理的軟骨細胞中分解代謝因子和軟骨基質(zhì)改變的調(diào)節(jié)中起關(guān)鍵作用〔27〕。趨化因子包含免疫細胞分泌的一小部分細胞因子或蛋白，且在免疫反應(yīng)過程中調(diào)節(jié)和控制免疫系統(tǒng)細胞〔28〕。鑒于趨化因子在炎癥期間白細胞運輸和活化過程中起關(guān)鍵性作用，趨化因子很可能在OA的發(fā)育發(fā)揮重要作用〔29〕。研究表明，在OA液中檢測到許多趨化因子可以促進細胞類型的轉(zhuǎn)化，包括CCL1848和CX3CL1〔30，31〕。