骨性關(guān)節(jié)炎滑膜核心基因篩選和生物信息學分析

張美齡 徐亞東 江小華

骨性關(guān)節(jié)炎(osteoarthritis, OA)是一種累及全關(guān)節(jié)的致殘性疾病，隨著全球人口老齡化和肥胖增加等因素的影響，發(fā)生率還在逐年上升,給個人健康及社會經(jīng)濟造成很大負擔[1]。OA主要病理特征包括關(guān)節(jié)軟骨退化、滑膜炎及關(guān)節(jié)結(jié)構(gòu)的改變[2]。目前OA的發(fā)病機制還不是完全清晰，關(guān)節(jié)置換是疾病終末期的有效治療方法，然而預后依然較差,因此疾病預防和早期OA診斷更為重要。此外，滑膜的病理改變可能發(fā)生在軟骨變化之前，滑膜組織在基因?qū)用娴淖兓赡苁羌膊“l(fā)生、發(fā)展的關(guān)鍵所在，隨著高通量微陣列技術(shù)的發(fā)展為深化基因-疾病關(guān)系的研究帶來了希望[3，4]。本研究通過對公共數(shù)據(jù)庫GEO來源的OA滑膜組織微陣列數(shù)據(jù)分析，包括篩選差異表達基因、功能富集分析、ROC曲線驗證，以期探索OA潛在診斷和治療靶點，同時為疾病分子機制研究提供生物信息學依據(jù)。

材料與方法

1.數(shù)據(jù)來源：從國際公開數(shù)據(jù)庫GEO (https:∥www.ncbi.nlm.nih.gov/geo/)下載骨性關(guān)節(jié)炎相關(guān)數(shù)據(jù)集GSE55235、GSE55457、GSE55584，物種選擇智人，3個微陣列數(shù)據(jù)集都基于GPL96平臺[5，6]。GSE55235包含10個OA滑膜組織樣本和10個正常滑膜組織樣本，GSE55457包含10個OA滑膜組織樣本和10個正常滑膜組織樣本，GSE55584包含6個OA滑膜組織樣本。

2.差異表達基因篩選：整合GSE55235、GSE55457、GSE55584微陣列數(shù)據(jù)集共26個OA滑膜組織樣本，20個正常滑膜組織樣本，并使用在線分析工具Network Analyst (https:∥www.networkanalyst.ca/)進行標準化處理，并進行主成分分析，再以P<0.05，|Log2(fold change, FC)|≥1為條件鑒定差異表達基因[7]。

3.差異表達基因的GO和KEGG通路分析：基因本體論(gene ontology, GO)術(shù)語包括生物過程(biological process, BP)、細胞組分(cellular components, CC)、分子功能(molecular function, MF)。使用R軟件(3.6.3版本)Org.Hs.eg.db包(ID轉(zhuǎn)換)、Ggplot2包(用于可視化)、Cluster-Profiler包(數(shù)據(jù)分析)進行GO和KEGG信號通路富集分析，調(diào)整后P<0.1為篩選標準。

4.蛋白-蛋白相互作用網(wǎng)絡(luò)構(gòu)建及核心蛋白篩選：使用在線分析工具STRING (https:∥string-db.org/)構(gòu)建基因編碼蛋白質(zhì)的相互作用(protein-protein interaction, PPI)網(wǎng)絡(luò)(可信度=0.4)并下載其TSV格式文件，標注基因上下調(diào)關(guān)系，利用細胞全景觀軟件Cyto-scape(3.8.3版本)展示DEGs的調(diào)控網(wǎng)絡(luò)，并使用軟件插件“Cytohubba”篩選最大相關(guān)性算法(maximal clique centrality，MCC)中位列前5個基因作為核心基因[8,9]。

5.核心基因的表達差異：使用在線分析工具Network Analyst進行核心基因的差異表達分析，以箱式圖形式展示。

6.核心基因診斷效能評價:使用R軟件(3.6.3版本)的PROC包(用于分析)、Ggplot2包(用于可視化)繪制5個核心基因的受試者工作特征(receiver operating characteristic curve, ROC)曲線。

結(jié) 果

1.數(shù)據(jù)預處理和差異表達基因篩選：整合3個芯片數(shù)據(jù)集共得到12545個基因的mRNA表達矩陣，并進行主成分分析顯示,OA組和對照組間樣本自然分開(圖1)。根據(jù)P<0.05，|Log2FC|≥1標準共鑒定出327個DEGs，包括219個上調(diào)基因和108個下調(diào)基因(圖2)。

圖1 微陣列數(shù)據(jù)集的主成分分析括號內(nèi)百分數(shù)表示OA組與對照組間差異中可以解釋全面分析結(jié)果的百分數(shù)，樣本與樣本之間的距離越小越相似，反之差異越大

圖2 差異表達基因火山圖每個點代表1個基因，紅色219個上調(diào)基因，藍色108個下調(diào)基因，灰色表示表達量差異無統(tǒng)計學意義的基因

2.差異表達基因的功能富集分析：通過Cluster-Profiler包富集分析，219個上調(diào)DEGs主要參與含膠原的細胞外基質(zhì)、質(zhì)膜外側(cè)、內(nèi)質(zhì)網(wǎng)腔的組成；參與細胞外基質(zhì)組織、白細胞遷移等生物過程；發(fā)揮信號蛋白受體激活、糖胺聚糖結(jié)合、胞外基質(zhì)結(jié)構(gòu)組成等分子功能，參與破骨細胞分化、軸突引導等KEGG通路(圖3A)。108個下調(diào)DEGs主要參與轉(zhuǎn)錄因子復合物的組成；參與糖皮質(zhì)激素反應、皮質(zhì)類固醇反應、脂多糖反應等生物過程；發(fā)揮DNA結(jié)合轉(zhuǎn)錄激活子激活、細胞因子激活、趨化因子激活等分子功能;富集于白細胞介素-17(interleukin-17, IL-17)、腫瘤壞死因子(tumor necrosis factor, TNF)、卡波西肉瘤相關(guān)皰疹病毒感染等KEGG通路(圖3B)。

圖3 差異表達基因GO和KEGG富集柱狀圖A.上調(diào)基因功能富集分析結(jié)果；B.下調(diào)基因功能富集分析結(jié)果

3.PPI網(wǎng)絡(luò)構(gòu)建及關(guān)鍵基因鑒定：基于在線分析工具STRING構(gòu)建(303個節(jié)點、873條邊組成)蛋白質(zhì)調(diào)控網(wǎng)絡(luò)，由Cyto-scape軟件進行PPI可視化展示(圖4A),利用Cytohubba插件MCC算法鑒定前5個核心基因分別是白介素-6(interleukin-6, IL-6)、C-X-C基序趨化因子配體8(c-x-c motif chemokine ligand8, CXCL8)、血管內(nèi)皮生長因子A(vascular endothelial growth factor A, VEGFA)、toll樣受體3(toll-like receptor, TLR3)、toll樣受體7(toll-like receptor, TLR7,圖4B)。

圖4 蛋白質(zhì)-蛋白質(zhì)相互作用網(wǎng)絡(luò)A.紅色表示上調(diào)的基因編碼蛋白質(zhì)，綠色表示下調(diào)的基因編碼蛋白質(zhì)；B.Cytohubba篩選5個核心基因編碼蛋白質(zhì)

4.核心基因的表達：使用Network Analyst數(shù)據(jù)庫分析發(fā)現(xiàn)，與對照組比較，TLR3和TLR7 滑膜mRNA水平在OA組表達升高(圖5A、B)；CXCL8、IL-6、VEGFA 滑膜mRNA水平則在OA組表達下降(圖5C～E)，以P<0.05為差異有統(tǒng)計學意義。

圖5 核心基因表達水平箱式圖A.TLR7；B.TLR3；C.CXCL8；D.IL-6；E.VEGFA

5.診斷效能評估：基于R語言繪制核心基因的ROC曲線其中IL-6(AUC=0.810)、CXCL8(AUC=0.846)、VEGFA(AUC=0.921)、TLR3基因(AUC=0.785)、TLR7基因(AUC=0.875，圖6)，AUC均>0.7具有較好的OA診斷效能。

圖6 核心基因診斷OA的ROC曲線A.表達上調(diào)基因TLR7、TLR3；B.表達下調(diào)基因CXCL8、IL-6、VEGFA

討 論

骨性關(guān)節(jié)炎是一種慢性退行性關(guān)節(jié)疾病，由機體體質(zhì)和機械等多因素復雜作用所致，常需要體格檢查、實驗結(jié)果、影像學來支持OA診斷[10，11]。基質(zhì)金屬蛋白酶(matrix metalloproteinases, MMPs)、C反應蛋白(C-reactive protein, CRP)、ⅡA型膠原氨基末端前肽(procollagen ⅡA amino-terminal pro-peptide, PⅡANP)等為目前報道有限數(shù)量的生物學標志物[12]。尚需探索更多高敏感度和特異性的診斷靶點來協(xié)助OA預防和治療，而且有證據(jù)表明滑膜炎發(fā)生在影像學變化之前，可能是OA發(fā)病和結(jié)構(gòu)進展的獨立驅(qū)動因素[3]。

本研究首先對3個GEO來源的骨性關(guān)節(jié)炎滑膜微陣列數(shù)據(jù)集進行主成分分析，發(fā)現(xiàn)骨性關(guān)節(jié)炎組和正常對照組間樣本自然分開，提示樣本來源可靠，接著構(gòu)建PPI網(wǎng)絡(luò)篩選得到327個差異表達基因，功能聚類分析發(fā)現(xiàn)DEGs主要集中于含膠原的細胞外基質(zhì)、白細胞遷移、細胞因子激活、趨化因子激活等基因本體術(shù)語；主要參與破骨細胞分化，軸突引導、IL-17、TNF等KEGG信號通絡(luò)。有實驗研究表明，在骨性關(guān)節(jié)炎初期，軟骨失去完整性，膠原外基質(zhì)降解，炎性介質(zhì)產(chǎn)生，并作用于臨近滑膜釋放IL-7、TNF-α等細胞因子[13,14]。這些結(jié)果與本研究發(fā)現(xiàn)具有較高吻合度，提示使用生物信息學方法可以從海量數(shù)據(jù)信息中，挖掘并尋找規(guī)律，為疾病研究提供一種較為高效的分析手段[15]。接著篩選關(guān)鍵調(diào)控基因以及ROC曲線驗證發(fā)現(xiàn)5個核心基因，其中TLR7、TLR3表達量上調(diào)，VEGFA、IL-6、CXCL8表達量下調(diào)，且AUC均>0.7(0.765～0.921)被認為具有較高診斷效能。

VEGFA為血管內(nèi)皮生長因子A,屬于VEGF家族中的一員，編碼肝素結(jié)合蛋白，對生理性和病理性血管生成均至關(guān)重要，VEGF以多種亞型形式存在，依賴于mRNA選擇性剪切。有研究顯示，VEGFA參與了OA病理過程，在有癥狀OA滑膜中基因和蛋白水平表達降低，具有促進血管生成和疼痛的作用[16,17]。IL-6作為廣泛表達于非特異性炎性介質(zhì)，在OA免疫、炎癥、軟骨破壞、軟骨重建中發(fā)揮重要作用[18]。有研究表明，IL-6在OA中似乎起著雙向調(diào)控作用，既能促進白細胞介素-1(interleukin-1, IL-1)、TNF受體、MMP阻滯劑下調(diào)炎性反應，又能提高免疫細胞功能和炎性反應，IL-6差異性表達至少說明OA存在慢性改變狀態(tài)[19]。CXCL8為C-X-C基序趨化因子配體8，該基因編碼的蛋白質(zhì)屬于CXC趨化家族的一員，是炎性反應的主要介質(zhì)，一項Meta分析顯示，CXCL8作為滑膜液生物學標志物具有較高的診斷價值[20]。本研究發(fā)現(xiàn)CXCL8 mRNA水平在OA組中表達下降，提示CXCL8表達量可能與OA嚴重程度相關(guān)，在疾病中可能潛在具有炎癥和穩(wěn)態(tài)雙重功能，其具體機制有待于進一步研究。

TLR3、TLR7編碼的蛋白是Toll樣受體家族成員，TLRs在病原體識別和固有免疫激活中發(fā)揮基本作用，其中TLR3可直接綁定β干擾素TIR結(jié)構(gòu)域銜接蛋白激活干擾素調(diào)節(jié)因子(interferon regulatory factor，IRFs)，進而產(chǎn)生炎性細胞因子、趨化因子和蛋白酶以促進OA病程進展。此外有研究表明OA與2個TLR3單核苷酸多態(tài)性(single nucleotide polymorphism, SNP)位點顯著相關(guān)，并且與TLR7的SNP位點也具有一定相關(guān)性，這些與OA風險遺傳相關(guān)的證據(jù)表明TLRs在OA發(fā)生、發(fā)展中具有重要作用。

綜上所述，本研究通過整合(GSE55235、GSE55584、GSE554573)微陣列數(shù)據(jù)集進行全基因組學分析，發(fā)現(xiàn)TLR7、TLR3、CXCL8、IL-6、VEGFA可能是骨性關(guān)節(jié)炎滑膜損害的關(guān)鍵基因，并初步探討了基因在疾病中的表達規(guī)律及診斷效能，這些都有助于闡明OA發(fā)生的分子機制，進一步基因功能聚類分析提示破骨細胞分化、細胞因子激活、TNF信號通路等可能參與了OA的發(fā)展過程。