基于生物信息與機器學習分析類風濕關節炎疾病特征基因與免疫細胞的關系①

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

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

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

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

1 資料與方法

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

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

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

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

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

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

2 結果

2.1 差異基因篩選 在GEO 網站［Home-GEONCBI（nih.gov）］下載原始數據，運用Perl 腳本對各基因進行標準化處理，處理后的數據用R 軟件進行再次處理，將GSE55235 和GSE12021 整理后的數據進行處理、合并、分析，共得到差異基因12 548 個，包含5 935 個上調基因和6 613 個下調基因，設定adj.P<0.05 和|log2FC|>2 為過濾條件，最終篩選出90 個差異基因，包含64 個上調基因和26 個下調基因，繪制差異基因熱圖和火山圖（圖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 圈圖。結果顯示GO 分析共得到主要條目209 個，KEGG 富集分析共得到11 個項目，其中GO BP 190 個，CC 5 個，MF 14 個，主要涉及白細胞激活、淋巴細胞活化、B 細胞受體信號通路、吞噬細胞介導的免疫、淋巴細胞分化、體液免疫反應等（圖2A）；KEGG 富集分析發現RA 相關通路有趨化因子信號通路、IL-17信號通路、Toll 樣受體信號通路、PPAR 信號通路、原發性免疫缺陷等（圖2B）。

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

2.3 機械學習法獲取疾病特征基因 運用LASSO回歸，通過R 軟件“glmnet”包對差異基因表達進行過濾，取誤差最小值的前10 個基因為疾病特征基因，繪制相關可視化圖（圖3A），縱坐標代表誤差大小，橫坐標代表基因個數，基因誤差按從高到低排名 依 次 得 到IGHM、SLAMF8、CXCL10、FNDC4、AIM2、EGR1、AKR1B10、FKBP5、TLR8、LIF，發 現第10個基因LIF誤差最小；運用SVM-RFE 法篩選疾病特征基因：通過R 軟件“e1071”“kernlab”“caret”包對差異基因表達進行過濾，交叉驗證尋找誤差最小的疾病特征基因并進行可視化分析（圖3B），縱坐標代表誤差大小，橫坐標代表基因個數，發現誤差表達最小的點共90個基因，將兩者各自得到的疾病特征基因取交集，繪制韋恩圖（圖4），得到IGHM、SLAMF8、CXCL10、FNDC4、AIM2、EGR1、AKR1B10、FKBP5、TLR8、LIF 10 個疾病特征基因，繪制10 個疾病特征基因的ROC 曲線，發現10 個疾病特征基因AUC 均>80%，運用測試基因芯片對上述10 個特征基 因 進 行 檢 測，發 現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 疾病特征基因與免疫細胞的相關性分析 使用 R 軟件（4.1.2 版本）“limma”“reshape2”“ggpubr”“ggExtra”包對5個疾病特征基因的免疫細胞浸潤結果與校正后的差異基因表達數據進行依次相關性分析，觀察各疾病特征基因表達與免疫細胞的相關性，繪制可視化棒棒糖圖（圖6），從左到右依次對應SLAMF8、IGHM、CXCL10、AIM2、AKR1B10，最左側縱坐標代表免疫細胞，最右側為P值大小，圓圈越大代表P絕對值越小，發現CD8 T 細胞、漿細胞、濾泡輔助細胞、激活記憶T 細胞CD4 在SLAMF8、IGHM、AIM2、CXCL10基因中的表達較高，且該基因與這些免疫細胞存在顯著相關性，活化的肥大細胞、記憶靜息CD4 T 細胞、單核細胞、激活的NK 細胞在AKR1B10基因中表達較高，且該基因與這些免疫細胞存在顯著相關性。SLAMF8還有M1巨噬細胞和B記憶細胞，IGHM 基因還有M1 巨噬細胞、M0 巨噬細胞和γδT 細胞，AIM2 還有M1 巨噬細胞與γδT 細胞，CXCL10還有M1巨噬細胞。

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

3 討論

RA 是一種慢性、多系統的自身免疫病，其特征為持續性滑膜炎癥，通常對稱分布于關節周圍，某些病例中出現全身癥狀。遺傳、環境因素和自身免疫3個重要因素在RA發病機制中發揮作用。

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

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

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

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