基于生物信息學(xué)探討骨關(guān)節(jié)炎相關(guān)分子機(jī)制及免疫細(xì)胞浸潤(rùn)分析

劉暢，卞華，許博

（1.河南省中醫(yī)院（河南中醫(yī)藥大學(xué)第二附屬醫(yī)院）風(fēng)濕病科，鄭州 450002；2.南陽(yáng)理工學(xué)院張仲景國(guó)醫(yī)國(guó)藥學(xué)院，河南省張仲景方藥與免疫調(diào)節(jié)重點(diǎn)實(shí)驗(yàn)室，河南 南陽(yáng) 473004）

骨關(guān)節(jié)炎（osteoarthritis，OA）可影響任何滑膜關(guān)節(jié)，其中髖部、膝部、手部、足部和脊柱是最常受累的部位［1］。目前對(duì)影響OA發(fā)病的危險(xiǎn)因素尚不清楚，多數(shù)研究［2］認(rèn)為可能與年齡、性別、肥胖、遺傳、飲食，以及關(guān)節(jié)損傷、錯(cuò)位和異常負(fù)荷有關(guān)。本研究基于基因表達(dá)圖譜數(shù)據(jù)庫(kù)芯片數(shù)據(jù)，通過(guò)與正常人群數(shù)據(jù)進(jìn)行對(duì)比，探討OA發(fā)病的分子機(jī)制。

1 材料與方法

1.1 數(shù)據(jù)收集

以“osteoarthritis”為關(guān)鍵詞，在基因表達(dá)綜合（gene expression omnibus，GEO）數(shù)據(jù)庫(kù)（https：//www.ncbi.nlm.nih.gov/geo/）中檢索，選取編號(hào)為GSE51588和GSE98918的2組芯片，共包含52例OA和22例正常對(duì)照組實(shí)驗(yàn)數(shù)據(jù)。再根據(jù)GEO數(shù)據(jù)庫(kù)中獲得的探針注釋文件，將每個(gè)數(shù)據(jù)集中的探針更改為基因符號(hào)，將這2個(gè)數(shù)據(jù)集合并到1個(gè)元數(shù)據(jù)隊(duì)列中進(jìn)行集成分析。當(dāng)同一基因符號(hào)對(duì)應(yīng)多個(gè)探針時(shí)，將探針的平均值作為該基因的表達(dá)值。然后選取包含12例OA和12例正常對(duì)照組樣本的GSE117999芯片作為驗(yàn)證隊(duì)列。

1.2 差異表達(dá)基因（differentially expressed gene，DEG）分析

使用R軟件中的“l(fā)imma”包對(duì)元數(shù)據(jù)隊(duì)列進(jìn)行差異分析，采用“FDR”處理方法并設(shè)置過(guò)濾條件為logFCfilter>1，校正P<0.05，進(jìn)行篩選得到DEG。然后使用R軟件中“clusterProfiler”包以P（Pvalue filter）<0.05為過(guò)濾條件對(duì)DEG進(jìn)行基因本體（Gene Ontology，GO）和疾病本體（Disease Ontology，DO）富集分析。

1.3 特征生物標(biāo)志物篩選

使用最小絕對(duì)收縮和選擇算子（least absolute shrinkage and selection operator，LASSO）與支持向量機(jī)（support vector machine，SVM）2種機(jī)器學(xué)習(xí)算法預(yù)測(cè)特征生物標(biāo)志物。首先利用R軟件中的“glmnet”包，對(duì)DEG使用LASSO回歸算法進(jìn)行交叉驗(yàn)證并篩選特征基因。通過(guò)R軟件中的“e1071”“kernlab”“caret”包采用支持向量機(jī)-遞歸特征消除（support vector machine-recursive feature elimination，SVM-RFE）方法，識(shí)別DEG中具有最高分辨能力的基因集。

1.4 特征生物標(biāo)志物對(duì)OA的診斷與治療價(jià)值判斷

為測(cè)試1.3中得到的特征基因價(jià)值，將驗(yàn)證組中所包含的12例OA和12例正常對(duì)照組樣本數(shù)據(jù)導(dǎo)入R軟件中，利用“l(fā)imma”和“ggpubr”包繪制箱線圖，判斷LASSO與SVM-RFE算法所得交集基因的價(jià)值。通過(guò)R軟件中的“pROC”包，使用相同的數(shù)據(jù)繪制受試者操作特征（receiver operating characteristic，ROC）曲線，用ROC曲線下面積（area under the curve，AUC）值確定OA的診斷有效性。

1.5 免疫細(xì)胞含量分析

采用CIBERSORT算法對(duì)浸潤(rùn)免疫細(xì)胞進(jìn)行相關(guān)分析和可視化。用R軟件中的“corrplot”和“vioplot”包對(duì)DEG行進(jìn)一步分析及可視化處理，得到OA免疫細(xì)胞之間的相關(guān)圖以及代表差異大小的小提琴圖。

1.6 已鑒定基因生物標(biāo)志物與浸潤(rùn)性免疫細(xì)胞的相關(guān)性分析

利用R軟件中的Spearman等級(jí)相關(guān)分析，研究已鑒定的基因生物標(biāo)志物與免疫細(xì)胞滲入水平的相關(guān)性，再用圖表函數(shù)和R軟件“ggplot2”包對(duì)二者的相關(guān)性進(jìn)行可視化處理。

2 結(jié)果

2.1 DEG獲取結(jié)果

對(duì)GSE51588、GSE98918芯片中基因表達(dá)數(shù)據(jù)進(jìn)行分析，結(jié)果顯示，共獲取OA組與正常對(duì)照組DEG 96個(gè)，其中上調(diào)基因47個(gè)，下調(diào)基因49個(gè)，包含顯著上調(diào)基因20個(gè)，顯著下調(diào)基因23個(gè)。見(jiàn)圖1。

圖1 訓(xùn)練組差異表達(dá)火山圖Fig.1 The training group showed differentially expressed volcano map

2.2 GO富集分析及DO富集分析結(jié)果

對(duì)DEG進(jìn)行GO和 DO富集分析，并繪制柱狀圖。GO富集分析結(jié)果顯示，DEG主要集中在共生體相互作用、真菌等生物學(xué)進(jìn)程和初級(jí)溶酶體、嗜苯胺藍(lán)顆粒等細(xì)胞組分，以及肝素結(jié)合、受體配體活動(dòng)等分子功能。DO富集結(jié)果顯示，DEG多富集在動(dòng)脈粥樣硬化、動(dòng)脈硬化性心血管疾病、動(dòng)脈硬化等疾病。

2.3 診斷特征生物標(biāo)志物的識(shí)別和驗(yàn)證

首先用LASSO算法縮小DEG的范圍，最終確定了19個(gè)變量作為OA的診斷生物標(biāo)志物（圖2A）。使用SVM-RFE算法確定了DEG中25個(gè)特征的子集（圖2B）。最終選擇了這2種算法之間13個(gè)重疊特征基因（圖2C）。

圖2 OA特征生物標(biāo)志物的識(shí)別Fig.2 Identification of OA characteristic biomarkers

為了產(chǎn)生更準(zhǔn)確可靠的結(jié)果，使用驗(yàn)證組GSE117999芯片驗(yàn)證這25個(gè)特征基因的表達(dá)水平，結(jié)果如圖3所示，樣品中CSN1S1、CXCL14、MTHFD2、NMNAT2、TLR7表達(dá)水平在OA中差異顯著（均P<0.05），因此，這5個(gè)顯著差異基因被用于logistic回歸算法建立診斷模型。

圖3 OA特征生物標(biāo)志物的驗(yàn)證Fig.3 Validation of OA characteristic biomarkers

2.4 特征基因的診斷有效性

如圖4所示，CSN1S1、CXCL14、MTHFD2、NMNAT2、TXR7的AUC值分別為0.819（95%CI：0.721～0.906）、0.899（95%CI：0.817～0.962）、0.862（95%CI：0.775～0.935）、0.884（95%CI：0.800～0.952）、0.858（95%CI：0.751～0.940），提示這5種特征基因生物標(biāo)志物對(duì)區(qū)分OA與正常對(duì)照組有統(tǒng)計(jì)學(xué)意義。

圖4 OA特征基因的診斷有效性Fig.4 Diagnostic validity of OA characteristic gene

2.5 免疫細(xì)胞浸潤(rùn)分析

OA組免疫細(xì)胞間相關(guān)性分析結(jié)果如圖5所示，正相關(guān)性較高的組合有靜息記憶CD4+T細(xì)胞和幼稚B細(xì)胞、幼稚B細(xì)胞和單核細(xì)胞、嗜酸性粒細(xì)胞和激活的樹(shù)突狀細(xì)胞等。負(fù)相關(guān)性較高的組合有幼稚CD4+T細(xì)胞和單核細(xì)胞、MO巨噬細(xì)胞和CD8+T細(xì)胞、中性粒細(xì)胞和幼稚CD4+T細(xì)胞等。免疫細(xì)胞差異分析結(jié)果如圖6所示，正常對(duì)照組中，幼稚B細(xì)胞、單核細(xì)胞、激活的肥大細(xì)胞、中性粒細(xì)胞的表達(dá)明顯增加（P<0.05）；OA組中，幼稚CD4+T細(xì)胞、濾泡輔助細(xì)胞、巨噬細(xì)胞M1、靜息樹(shù)突狀細(xì)胞的表達(dá)明顯增加（P<0.05）。

圖5 免疫細(xì)胞間的相關(guān)性Fig.5 The correlation between immune cells

圖6 免疫細(xì)胞差異分析Fig.6 Analysis of immune cell differences

2.6 基因生物標(biāo)志物與浸潤(rùn)性免疫細(xì)胞的相關(guān)性分析

相關(guān)性分析結(jié)果顯示，CSN1S1與單核細(xì)胞、漿細(xì)胞呈負(fù)相關(guān)（圖7A）；MTHFD2與靜息肥大細(xì)胞呈正相關(guān)，與CD8+T細(xì)胞呈負(fù)相關(guān)（圖7B）；NMNAT2與巨噬細(xì)胞M1、調(diào)節(jié)性T細(xì)胞呈正相關(guān)，與幼稚B細(xì)胞呈負(fù)相關(guān)（圖7C）。

圖7 基因生物標(biāo)志物與浸潤(rùn)性免疫細(xì)胞的相關(guān)性Fig.7 Association of genetic biomarkers with invasive immune cells

3 討論

OA是一種慢性退行性疾病，常見(jiàn)于40歲以上人群［3］，是致殘的主要原因［4］。OA的病理特征是滑膜關(guān)節(jié)中關(guān)節(jié)軟骨的局灶性缺失，與不同程度的骨贅形成、軟骨下骨改變和滑膜炎相關(guān)［5］。雖然滑膜炎癥與疼痛嚴(yán)重程度密切相關(guān)，但導(dǎo)致OA疼痛的分子機(jī)制尚不清楚。本研究通過(guò)生物信息學(xué)技術(shù)對(duì)比OA組與正常對(duì)照組后，共篩查出96個(gè)DEG。GO和DO富集分析顯示，這些DEG涉及多種生物學(xué)進(jìn)程、細(xì)胞組分、分子功能及疾病。

本研究中，基于LASSO與SVM這2種機(jī)器學(xué)習(xí)算法，還識(shí)別出CSN1S1、CXCL14、MTHFD2、NMNAT2、TLR7共5種特征基因。BROPHY等［6］通過(guò)生物信息學(xué)技術(shù)在OA和非OA樣本之間鑒定出168個(gè)顯著差異表達(dá)轉(zhuǎn)錄物，發(fā)現(xiàn)CSN1S1是OA半月板中升高較突出的轉(zhuǎn)錄物。LIU等［7］用激光捕獲顯微解剖1周齡小鼠脛骨關(guān)節(jié)軟骨生長(zhǎng)的淺層、中層和深層區(qū)域并分離軟骨細(xì)胞，通過(guò)微陣列基因表達(dá)模式鑒定出中層關(guān)節(jié)軟骨中含有CXCL14等分子標(biāo)志物。ZHAO等［8］通過(guò)差異基因表達(dá)分析篩選出MTHFD2、TLR7、CX3CR1、SLC7A5、ARL4C 5種生物標(biāo)志物，在軟骨、軟骨下骨和滑膜組織中表達(dá)差異顯著，并證明使用這5種生物標(biāo)志物通過(guò)邏輯回歸合成的模型診斷OA的準(zhǔn)確性優(yōu)于單一生物標(biāo)志物，可作為新型生物標(biāo)志物提高膝OA診斷的準(zhǔn)確度，改善臨床療效。HOSHIKAWA等［9］通過(guò)OA大鼠模型研究發(fā)現(xiàn)，滑膜組織釋放的細(xì)胞外miR-21通過(guò)激活TLR7介導(dǎo)OA模型大鼠的膝關(guān)節(jié)疼痛，miR-21的鎮(zhèn)痛作用可通過(guò)拮抗TLR7被阻斷，TLR7拮抗劑對(duì)膝OA疼痛有持久的鎮(zhèn)痛作用。TLR7參與免疫反應(yīng)已在多種炎癥性疾病中報(bào)道，OA患者的軟骨組織和脂多糖誘導(dǎo)的軟骨細(xì)胞中TLR7表達(dá)增加，TLR7下調(diào)可通過(guò)阻斷p21介導(dǎo)的JAK2/STAT3途徑，減弱軟骨細(xì)胞的凋亡和損傷，表明TLR7可能是OA的一個(gè)治療靶點(diǎn)［10］。研究［11］發(fā)現(xiàn)，TLR7通過(guò)不同免疫細(xì)胞、多種途徑影響OA免疫微環(huán)境，并在OA患者中表達(dá)水平較高，可作為OA的診斷標(biāo)志物。另有研究［12］通過(guò)SNaPshot測(cè)序技術(shù)證明，TLR7可能在中國(guó)漢族人群的膝OA中發(fā)揮易感性作用，并可能與膝OA的嚴(yán)重程度和膝OA積液滑膜炎的風(fēng)險(xiǎn)有關(guān)。

免疫細(xì)胞浸潤(rùn)分析結(jié)果顯示，特征基因與單核細(xì)胞、巨噬細(xì)胞M1、漿細(xì)胞、CD8+T細(xì)胞、幼稚B細(xì)胞、調(diào)節(jié)性T細(xì)胞、靜息肥大細(xì)胞等相關(guān)。單核細(xì)胞、巨噬細(xì)胞是健康關(guān)節(jié)的主要免疫細(xì)胞［13］，參與調(diào)節(jié)骨轉(zhuǎn)換，清除關(guān)節(jié)腔中的碎片，并與滑膜成纖維細(xì)胞一起形成保護(hù)屏障，是創(chuàng)傷后OA的治療靶點(diǎn)。無(wú)論是巨噬細(xì)胞M1亞型還巨噬細(xì)胞M2亞型，均與OA的發(fā)病機(jī)制有關(guān)［14］。OA小鼠模型滑膜巨噬細(xì)胞M1極化部分通過(guò)分泌 Rspo2 和激活軟骨細(xì)胞中的β-連環(huán)蛋白信號(hào)傳導(dǎo)，導(dǎo)致小鼠OA惡化，提示巨噬細(xì)胞M1和 Rspo2 是OA的潛在治療靶點(diǎn)［15-16］。果苷可通過(guò)抑制體外和體內(nèi)巨噬細(xì)胞M1極化，降低白細(xì)胞介素（interleukin，IL）-6 和 TNF-α的分泌，延緩軟骨和細(xì)胞外基質(zhì)降解和軟骨細(xì)胞肥大，從而減輕OA小鼠的滑膜炎癥并延緩OA的發(fā)展。研究［17］顯示，漿細(xì)胞和漿母細(xì)胞可產(chǎn)生如IL-17和IL-10等細(xì)胞因子，在自身免疫性疾病中發(fā)揮調(diào)節(jié)作用。另有研究［18］顯示，富含血小板的血漿由于其軟骨誘導(dǎo)特性和豐富的生長(zhǎng)因子庫(kù)而被廣泛用于治療軟骨缺損和OA，臨床上對(duì)OA患者有很明確的療效。在OA發(fā)病過(guò)程中，CD8+T細(xì)胞被激活并在OA進(jìn)展過(guò)程中擴(kuò)增，抑制關(guān)節(jié)中TIMP-1的表達(dá)，從而延緩OA病情進(jìn)展［19］。調(diào)節(jié)性T細(xì)胞通過(guò)促進(jìn)抗炎反應(yīng)之間的最佳平衡，產(chǎn)生免疫生理適應(yīng)，包括減少全身促炎癥介質(zhì)產(chǎn)生，刺激由中性粒細(xì)胞介導(dǎo)的針對(duì)病原體的先天防御，即吞噬作用，實(shí)現(xiàn)最佳的先天反應(yīng)［20］。肥大細(xì)胞可通過(guò)淋巴器官內(nèi)的同源相互作用微調(diào)T細(xì)胞和B細(xì)胞活性，從而調(diào)節(jié)局部及全身炎癥反應(yīng)［21］。

綜上所述，本研究發(fā)現(xiàn)，OA發(fā)病可能是通過(guò)調(diào)節(jié)CSN1S1、CXCL14、MTHFD2、NMNAT2、TLR7等特征基因表達(dá)，參與單核細(xì)胞、巨噬細(xì)胞M1、漿細(xì)胞、CD8+T細(xì)胞、幼稚B細(xì)胞、調(diào)節(jié)性T細(xì)胞、靜息肥大細(xì)胞等實(shí)現(xiàn)。本研究存在數(shù)據(jù)庫(kù)樣本數(shù)量少、未行細(xì)胞分子學(xué)實(shí)驗(yàn)驗(yàn)證等缺點(diǎn)，可能存在誤差。總之，本研究應(yīng)用生物信息學(xué)技術(shù)探討了OA發(fā)病的相關(guān)分子機(jī)制，為今后OA的治療提供了新的靶點(diǎn)視角。