宋世雷 陳躍平 陳 鋒 （廣西中醫(yī)藥大學附屬瑞康醫(yī)院，南寧 530011）

類風濕關(guān)節(jié)炎（rheumatoid arthritis，RA）是以關(guān)節(jié)滑膜慢性炎癥為主要病理特征的一種自身免疫性疾病，其發(fā)展與自身抗體病理過程相關(guān)，發(fā)病初期可在血清中發(fā)現(xiàn)類風濕因子（rheumatoid factor，RF）和抗瓜氨酸蛋白抗體（anti-citrulline protein antibody，ACPA）水平升高［1］。隨著各種免疫細胞與成纖維細胞和細胞因子相互作用，滑膜組織逐漸產(chǎn)生慢性炎癥，并發(fā)生骨侵蝕和關(guān)節(jié)破壞，出現(xiàn)臨床癥狀和損傷，進而發(fā)展為嚴重殘疾，此外全身炎癥可損害多個器官系統(tǒng)，如心臟組織、血管系統(tǒng)、腎臟、肺組織和神經(jīng)系統(tǒng)，嚴重影響患者生活質(zhì)量［2-3］。

RA 發(fā)病機制較為復雜，目前普遍認為其發(fā)病機制與遺傳、表觀遺傳學、環(huán)境、代謝、免疫和微生物因素相互作用有關(guān)，尚未得到系統(tǒng)性闡述。臨床研究和實驗室研究均發(fā)現(xiàn)RA 疾病進展受多種免疫細胞影響，如T 細胞、B 細胞、單核細胞、巨噬細胞、中性粒細胞、肥大細胞、樹突狀細胞、Treg 細胞和纖維細胞樣滑膜細胞等［4］。這些免疫細胞在疾病微環(huán)境中表現(xiàn)出不同作用，而作用不同取決于疾病環(huán)境與疾病特征基因表達不同，因此免疫細胞介導的滑膜炎癥和軟骨破壞機制可能在不同患者中有不同表現(xiàn)，導致RA 存在臨床異質(zhì)性，極大影響治療的準確性［5］。因此需要尋找一種可通過檢測不同患者疾病特征基因表達情況進行提前診斷與評估患者RA發(fā)生及發(fā)生后免疫細胞浸潤情況的方法，實現(xiàn)早期診斷與精準治療。

生物信息學作為一種新型研究手段，可高效系統(tǒng)地了解疾病分子作用方式，利用公共基因芯片數(shù)據(jù)庫中的RA 患者與健康者滑膜組織進行合并分析，尋找RA 滑膜組織中的差異基因與分子通路，運用機械學習法篩選RA 疾病特征基因并探討疾病特征基因與免疫浸潤情況的相關(guān)性，為系統(tǒng)闡述發(fā)病機制及預防和提早診斷RA 提供更多依據(jù)，最終實現(xiàn)精準治療。

1 資料與方法

1.1 數(shù)據(jù)來源 從GEO 數(shù)據(jù)庫下載GSE12021、GSE55235 和 GSE77298 基因表達譜芯片。GSE12021 數(shù)據(jù)集包含12 例RA 患者和9 例正常關(guān)節(jié)患者滑膜組織，RA 滑膜組織由德國漢諾威臨床化學研究所胡貝兒教授對該類患者行關(guān)節(jié)置換術(shù)或滑膜切除術(shù)后提取，正常滑膜組織則由創(chuàng)傷性關(guān)節(jié)骨折患者提?。?］。GSE55235數(shù)據(jù)集包括10例RA患者和10 例正常關(guān)節(jié)患者滑膜組織，RA 滑膜組織由WOETZEL 等［7］在進行關(guān)節(jié)置換術(shù)或滑膜切除術(shù)時收集，正常組滑膜組織來源于死亡患者正常關(guān)節(jié)或創(chuàng)傷性關(guān)節(jié)骨折患者；GSE77298數(shù)據(jù)集包括16例RA患者和7例正常關(guān)節(jié)患者滑膜組織，由Mathijs活檢 所 得，基 于GLP570 平 臺 完 成［8］；GSE12021、GSE55235數(shù)據(jù)集均基于GPL96平臺完成。

1.2 數(shù)據(jù)處理及差異基因篩選 運用Perl 腳本對GSE12021、GSE55235 數(shù)據(jù)集進行整理，R 軟件（版本4.1.2）“SVA”包對整理后的數(shù)據(jù)進行處理合并，獲取基因批次信息，運用R軟件“l(fā)imma”包分別對兩組合并后的批次信息進行常規(guī)矯正，輸出矯正后的結(jié)果，篩選矯正后同時符合adj.P<0.05 和|log2FC|>2 的基因作為差異表達基因，分析得到差異基因，輸出矯正后的差異基因表達和差異結(jié)果，“pheatmap”“ggplot2”軟件包構(gòu)建整合數(shù)據(jù)集中差異表達基因的分層聚類圖與火山圖。

1.3 差異表達基因GO和KEGG富集分析 為探求差異表達基因主要功能和途徑，通過R 軟件“Bioconductor”包對差異表達基因進行GO 和KEGG 富集分析，篩選有顯著差異的基因（P<0.05）；GO 富集分析包含生物學過程（biological process，BP）、細胞成分（cellular constituent，CC）和分子功能（molecular function，MF），通過“ggplot2”包對富集結(jié)果進行可視化，分別繪制GO和KEGG圈圖、氣泡圖和柱狀圖。

1.4 機器學習法獲取疾病特征基因 運用機器學習法在差異基因中尋找疾病特征基因，首先運用LASSO回歸法：R 軟件“glmnet”包對差異基因表達進行過濾，通過交叉驗證尋找誤差最小的基因則為疾病特征基因。運用SVM-RFE 法篩選疾病特征基因：R 軟件“e1071”“kernlab”“caret”包對差異基因表達進行過濾，通過交叉驗證尋找誤差最小的基因則為疾病特征基因，將兩者篩選的疾病特征基因取交集，繪制VENN圖，篩選最終的疾病特征基因，對篩選獲得的基因進行ROC 曲線檢測，GSE77298 與交集疾病特征基因進行驗證，繪制ROC曲線。

1.5 芯片的免疫浸潤分析 使用R 軟件（4.1.2 版本）通過CIBERSORT 算法對以上獲得的聯(lián)合芯片的矯正基因表達矩陣進行分析，篩選P<0.05 的樣品，并計算各樣品中22種免疫細胞占比。利用“vioplot”軟件包比較篩選的RA 樣品與正常樣品滑膜組織，分析兩組間各免疫細胞浸潤水平，P<0.05 代表兩組間差異有統(tǒng)計學意義。

1.6 疾病特征基因與免疫細胞的關(guān)系 使用R 軟件（4.1.2 版 本）“l(fā)imma”“reshape2”“ggpubr”“ggExtra”包對疾病特征基因免疫細胞浸潤結(jié)果與校正后的差異基因表達數(shù)據(jù)進行依次相關(guān)性分析，觀察各疾病特征基因表達與免疫細胞的相關(guān)性，并對其統(tǒng)計結(jié)果進行可視化分析，繪制相關(guān)性棒棒糖圖。

2 結(jié)果

2.1 差異基因篩選 在GEO 網(wǎng)站［Home-GEONCBI（nih.gov）］下載原始數(shù)據(jù)，運用Perl 腳本對各基因進行標準化處理，處理后的數(shù)據(jù)用R 軟件進行再次處理，將GSE55235 和GSE12021 整理后的數(shù)據(jù)進行處理、合并、分析，共得到差異基因12 548 個，包含5 935 個上調(diào)基因和6 613 個下調(diào)基因，設定adj.P<0.05 和|log2FC|>2 為過濾條件，最終篩選出90 個差異基因，包含64 個上調(diào)基因和26 個下調(diào)基因，繪制差異基因熱圖和火山圖（圖1）。圖1A 中正常組為Control 組，疾病組為Treat 組，紅色表示高表達，藍色為低表達，顏色越深代表程度越高，圖1B紅色為高表達，綠色為低表達。

圖1 差異基因熱圖與火山圖Fig.1 Heatmap and volcanic map of difference genes

2.2 差異基因GO與KEGG 富集分析 通過R軟件（4.1.2）對篩選出的90 個差異基因進行GO 和KEGG 富集分析，得到GO 和KEGG 圈圖。結(jié)果顯示GO 分析共得到主要條目209 個，KEGG 富集分析共得到11 個項目，其中GO BP 190 個，CC 5 個，MF 14 個，主要涉及白細胞激活、淋巴細胞活化、B 細胞受體信號通路、吞噬細胞介導的免疫、淋巴細胞分化、體液免疫反應等（圖2A）；KEGG 富集分析發(fā)現(xiàn)RA 相關(guān)通路有趨化因子信號通路、IL-17信號通路、Toll 樣受體信號通路、PPAR 信號通路、原發(fā)性免疫缺陷等（圖2B）。

圖2 GO與KEGG圈圖Fig.2 GO and KEGG circle graphs

2.3 機械學習法獲取疾病特征基因 運用LASSO回歸，通過R 軟件“glmnet”包對差異基因表達進行過濾，取誤差最小值的前10 個基因為疾病特征基因，繪制相關(guān)可視化圖（圖3A），縱坐標代表誤差大小，橫坐標代表基因個數(shù)，基因誤差按從高到低排名 依 次 得 到IGHM、SLAMF8、CXCL10、FNDC4、AIM2、EGR1、AKR1B10、FKBP5、TLR8、LIF，發(fā) 現(xiàn)第10個基因LIF誤差最??；運用SVM-RFE 法篩選疾病特征基因：通過R 軟件“e1071”“kernlab”“caret”包對差異基因表達進行過濾，交叉驗證尋找誤差最小的疾病特征基因并進行可視化分析（圖3B），縱坐標代表誤差大小，橫坐標代表基因個數(shù)，發(fā)現(xiàn)誤差表達最小的點共90個基因，將兩者各自得到的疾病特征基因取交集，繪制韋恩圖（圖4），得到IGHM、SLAMF8、CXCL10、FNDC4、AIM2、EGR1、AKR1B10、FKBP5、TLR8、LIF 10 個疾病特征基因，繪制10 個疾病特征基因的ROC 曲線，發(fā)現(xiàn)10 個疾病特征基因AUC 均>80%，運用測試基因芯片對上述10 個特征基 因 進 行 檢 測，發(fā) 現(xiàn)IGHM、SLAMF8、CXCL10、AIM2、AKR1B10 AUC均>0.8，說明IGHM、SLAMF8、CXCL10、AIM2、AKR1B10為疾病特征基因。

圖3 LASSO回歸圖與SVM-RFE圖Fig.3 LASSO regression and SVM-RFE

圖4 VENN圖Fig.4 VENN diagram

2.4 芯片免疫浸潤分析 CIBERSORT 算法分析矯正后的兩個芯片基因表達，其中RA 組有22 個滑膜組織樣本，正常組有19 個滑膜組織樣本，同時篩選得到符合條件的22 個RA 與19 個正常組的22 個免疫細胞浸潤存在差異，計算41個樣品免疫細胞占比（圖5A）。與正?；は啾?，RA 滑膜中漿細胞、CD8 T 細胞、記憶靜息T 細胞CD4、記憶激活T 細胞CD4、T細胞濾泡輔助器、激活的NK 細胞、單核細胞、巨噬細胞M1、活化的肥大細胞9 種免疫細胞與正常組織存在明顯差異（圖5B）。

圖5 免疫細胞占比柱狀圖及小提琴圖Fig.5 Proportion of immune cells columnar and vioplot

2.5 疾病特征基因與免疫細胞的相關(guān)性分析 使用 R 軟件（4.1.2 版本）“l(fā)imma”“reshape2”“ggpubr”“ggExtra”包對5個疾病特征基因的免疫細胞浸潤結(jié)果與校正后的差異基因表達數(shù)據(jù)進行依次相關(guān)性分析，觀察各疾病特征基因表達與免疫細胞的相關(guān)性，繪制可視化棒棒糖圖（圖6），從左到右依次對應SLAMF8、IGHM、CXCL10、AIM2、AKR1B10，最左側(cè)縱坐標代表免疫細胞，最右側(cè)為P值大小，圓圈越大代表P絕對值越小，發(fā)現(xiàn)CD8 T 細胞、漿細胞、濾泡輔助細胞、激活記憶T 細胞CD4 在SLAMF8、IGHM、AIM2、CXCL10基因中的表達較高，且該基因與這些免疫細胞存在顯著相關(guān)性，活化的肥大細胞、記憶靜息CD4 T 細胞、單核細胞、激活的NK 細胞在AKR1B10基因中表達較高，且該基因與這些免疫細胞存在顯著相關(guān)性。SLAMF8還有M1巨噬細胞和B記憶細胞，IGHM 基因還有M1 巨噬細胞、M0 巨噬細胞和γδT 細胞，AIM2 還有M1 巨噬細胞與γδT 細胞，CXCL10還有M1巨噬細胞。

圖6 疾病特征基因與免疫細胞相關(guān)性的棒棒糖圖Fig.6 Lollipop chart of correlation between disease characteristic genes and immune cells

3 討論

RA 是一種慢性、多系統(tǒng)的自身免疫病，其特征為持續(xù)性滑膜炎癥，通常對稱分布于關(guān)節(jié)周圍，某些病例中出現(xiàn)全身癥狀。遺傳、環(huán)境因素和自身免疫3個重要因素在RA發(fā)病機制中發(fā)揮作用。

通過公共基因芯片數(shù)據(jù)庫網(wǎng)站下載RA 芯片數(shù)據(jù)，對其進行數(shù)據(jù)合并分析，設定篩選條件后共篩選出差異基因90 個，包含64 個上調(diào)基因和26 個下調(diào)基因，對上述基因進行GO 與KEGG 富集分析發(fā)現(xiàn)，GO 結(jié)果顯示其主要富集于淋巴細胞活化、B 細胞受體信號通路、吞噬細胞介導的免疫、淋巴細胞分化、體液免疫反應，上述機體功能表達均與免疫息息相關(guān)［9］。淋巴細胞是白細胞中最小的一種免疫細胞，由淋巴器官產(chǎn)生，根據(jù)淋巴細胞來源、表面標志物與免疫功能主要分為T 淋巴細胞、B 淋巴細胞和NK 細胞，研究發(fā)現(xiàn)以上3 種免疫細胞在RA 發(fā)生機制上有各自的重要作用［10-12］。當T 淋巴細胞被激活，經(jīng)歷一系列轉(zhuǎn)變得到致敏T淋巴細胞，發(fā)揮抗感染、抗腫瘤、抗異體細胞作用，T 淋巴細胞激活和分化是促進炎癥發(fā)生發(fā)展的重要生理病理過程。研究發(fā)現(xiàn)RA 初始階段涉及T 細胞和B 細胞激活，由T細胞與B細胞主導和組織促炎細胞因子參與炎癥反應，這些炎癥因子包括TNF-α、IL-1、IL-17等，最終刺激骨和軟骨炎癥及降解［13］。B 淋巴細胞主要通過產(chǎn)生并分泌免疫球蛋白和細胞內(nèi)因子，包括抗原呈遞和T細胞激活，與趨化因子相互作用，促進TNF-α形成，激活巨噬細胞，參與免疫調(diào)節(jié)［14］。另有研究認為B細胞在RA發(fā)病機制中發(fā)揮更重要作用，主要原因為B 細胞產(chǎn)生RF（抗IgG 自身抗體），且B 細胞介導的抗CCP 產(chǎn)生后，RA 診斷的特異度已提高至90%～98%［15］。KEGG 富集分析發(fā)現(xiàn)RA 相關(guān)通路有趨化因子信號通路、IL-17 信號通路、Toll 樣受體信號通路、PPAR 信號通路、原發(fā)性免疫缺陷等。IL-17屬于T 細胞子集，主要由CD4+Th17 細胞產(chǎn)生，也可由CD8+T 細胞、NKT 細胞、γδT 細胞、中性粒細胞和淋巴組織誘導劑樣細胞產(chǎn)生，IL-17 家族由IL-17A～IL-17F 6 個成員組成，參與自身免疫病發(fā)生發(fā)展［16］。多個動物實驗發(fā)現(xiàn)IL-17 信號通路與RA 顯著相關(guān)，主要與炎癥發(fā)生、軟骨細胞破壞、骨質(zhì)破壞有關(guān)［17-19］。有證據(jù)表明IL-17 信號通路是關(guān)節(jié)炎癥發(fā)生的重要誘因，在RA 大鼠關(guān)節(jié)中注射IL-17 后發(fā)現(xiàn)中性粒細胞增殖與功能增強，可能與IL-17 信號通過誘導各種趨化因子和多種炎癥介質(zhì)分泌，包括IL-8、CXCL1、CXCL6、IL-1β、IL-6、TNF-α、GM-CSF、MIP-2、MCP-1 和G-CSF 等，促進中性粒細胞和單核細胞募集有關(guān)，IL-17 還發(fā)揮IL-17RA 和IL-17RC 介導的抗凋亡作用，促進滑膜組織增生，表明IL-17 通過刺激滑膜炎癥和滑膜組織增生引發(fā)RA［20］。IL-17還可上調(diào)RANK 和M-CSF 因子，在RANKL 和M-CSF刺激下，破骨細胞分泌各種中間體和促炎介質(zhì)，如誘導型一氧化氮合酶，促進骨丟失，可見IL-17 信號通路在RA 發(fā)生中占據(jù)重要位置［21］。Toll 樣受體屬于蛋白質(zhì)識別受體，人類中Toll 樣受體家族由10 個成員組成，即TLR1～TLR10，Toll 樣受體通過感知保守的病原體相關(guān)微生物模式和細胞損傷后產(chǎn)生的危險相關(guān)分子模式變化引發(fā)固有免疫應答［22］。研究發(fā)現(xiàn)RA 中TLR2、TLR3、TLR4、TLR7 免疫活性在滑膜組織的內(nèi)膜和下層表達上調(diào)，而在骨關(guān)節(jié)炎患者和健康供體滑膜組織中低表達，早期RA 滑膜中主要表現(xiàn)為TLR3 和TLR4 表達上調(diào)，而后期RA 滑膜則表現(xiàn)為TLR2、TLR3 和TLR4 上調(diào)，Toll 樣受體特異性高表達說明其在RA 發(fā)病中的重要作用［22］。Toll 樣受體信號通路對RA 的作用機制主要通過與外源性或內(nèi)源性配體結(jié)合被激活，其中髓系分化原代反應蛋白88（MyD88）和誘導含TLR 結(jié)構(gòu)域的適配器蛋白IFN-β（TRIF）是TLR 信號轉(zhuǎn)導的兩個主要配體，TLRs與MyD88結(jié)合后刺激IL-1受體相關(guān)激酶（IRAK）磷酸化，導致TRAF（TNF 受體相關(guān)因子）募集，從而激活NF-κB 信號通路和MAPK 信號通路，MAPK 信號通路下游包括細胞外信號調(diào)節(jié)激酶（ERK1/2）、C-Junn 端激酶（JNK）和P38，最終誘導促炎細胞因子TNF-α、IL-1 和IL-12 表達，提高炎癥細胞趨化因子與IFN-α/β 誘導的基因產(chǎn)物表達，TLRs與TRIF結(jié)合則誘導IRF3/7表達從而引起IFN-β表達升高，說明Toll樣受體信號通路是RA炎癥發(fā)生的機制之一［23-24］。RA 病理發(fā)展特征為滑膜腫瘤樣擴張，隨后鄰近關(guān)節(jié)軟骨和骨被破壞，主要原因為活化的纖維細胞樣滑膜細胞過度增殖、遷移和激活，侵襲鄰近軟骨和營養(yǎng)骨血管，因此抑制滑膜細胞增殖和遷移是治療RA 的理想靶點［25-27］。研究表明PPAR-γ可能參與滑膜細胞表達調(diào)控，RA 成纖維樣滑膜細胞中，PPAR-γ 表達下調(diào)，且PPAR-γ 抑制劑增加滑膜細胞增殖和遷移，而PPAR-γ激動劑表達抑制滑膜細胞增殖和遷移，為后續(xù)通過PPAR-γ 信號通路對RA 進行治療提供了新思路［28-29］。GO 與KEGG 富集分析發(fā)現(xiàn)RA 發(fā)生與淋巴細胞、IL-17 信號通路、Toll樣受體信號通路激活促進炎癥因子釋放及PPAR 信號通路調(diào)控纖維細胞樣滑膜細胞密切相關(guān)。

機器學習是一門多領(lǐng)域交叉學科，涉及概率論、統(tǒng)計學、逼近論、凸分析、算法復雜度理論等多學科，其常見算法包括決策樹算法、樸素貝葉斯算法、支持向量機算法、隨機森林算法和人工神經(jīng)網(wǎng)絡算法等［30］。本研究運用生物信息學中強大的特征選擇算法（SVM-RFE 和LASSO 回歸法）分別處理芯片數(shù)據(jù)，利用SVM 在數(shù)據(jù)集訓練得到的權(quán)重向量對特征進行排序，剔除無用特征，用剩余特征再次訓練模型進行下一次迭代，選出需要的特征數(shù)，保留交叉驗證后誤差最小的疾病特征基因，運用LASSO回歸法將處理的差異表達基因進行特征篩選，剔除無相關(guān)性或相關(guān)性不顯著的疾病特征基因，得到相關(guān)性顯著的疾病特征基因。將兩種算法得出的結(jié)果取交集，并對所得結(jié)果進行檢驗，繪制ROC 曲線，用檢測基因樣本對得到的結(jié)果進行檢驗，最終得到AUC>80%的5 個基因為疾病特征基因，即SLAMF8、IGHM、CXCL10、AIM2、AKR1B10。

免疫細胞浸潤是炎癥反應的前提，免疫細胞與RA 基因表達可能存在密切聯(lián)系，研究發(fā)現(xiàn)調(diào)節(jié)免疫細胞浸潤情況可有效控制RA 發(fā)生。因此本研究通過CIBERSORT 算法分析免疫細胞在RA 滑膜中的表達，并通過將篩選的疾病特征基因表達與免疫浸潤表達進行相關(guān)性分析，結(jié)果發(fā)現(xiàn)RA 滑膜中漿細胞、CD8 T 細胞、記憶靜息T 細胞CD4、記憶激活T細胞CD4、T 細胞濾泡輔助器、激活的NK 細胞、單核細胞、M1 巨噬細胞、活化的肥大細胞9 種免疫細胞與正?；ご嬖诿黠@差異，這種差異與疾病特征基因顯著相關(guān)，發(fā)現(xiàn)CD8 T 細胞、漿細胞、濾泡輔助細胞、激活記憶T 細胞CD4 在SLAMF8、IGHM、AIM2、CXCL10 基因中表達較高，且該基因與這些免疫細胞顯著相關(guān)，活化的肥大細胞、記憶靜息CD4 T 細胞、單核細胞、激活的NK細胞在AKR1B10基因中表達較高，且該基因與這些免疫細胞顯著相關(guān)。此外SLAMF8 還有M1 巨噬細胞和B 記憶細胞，IGHM 基因還有M1 巨噬細胞、M0 巨噬細胞和γδT 細胞，AIM2 還有M1 巨噬細胞與γδT 細胞，CXCL10 還有M1 巨噬細胞，不僅證實了5 個疾病特征基因?qū)A免疫應答的重要性，也可以看出不同基因影響不同免疫細胞在體內(nèi)的表達，從而影響RA發(fā)生。

綜上，本研究通過生信分析與機器學習法得到了RA 的主要疾病特征基因及免疫浸潤情況，并闡明了不同特征基因與免疫細胞表達的相關(guān)性，說明可通過判斷患者特征基因表達預測RA 免疫細胞表達，對不同患者進行針對性治療，提高疾病治愈率，提供臨床診療思路。