基于機器學習篩選類風濕關節炎的診斷標志基因和免疫浸潤分析

李玲琴, 周睿姣, 張燕妮, 賀泓霓, 袁心柱

類風濕關節炎(rheumatoid arthritis,RA)是一種慢性炎癥性自身免疫疾病[1],主要影響滑膜關節,表現為滑膜微血管生成、滑膜襯里細胞增生,滑膜間質有大量免疫細胞浸潤刺激局部炎癥,若不及時治療,可能會導致永久性關節損傷和殘疾[2]。全世界大約每200名成年人中就有1人受RA的影響,且女性的發病率是男性的2～3倍,發病高峰期為50～59歲。據統計,RA患者的壽命較無RA者少6～11年[3]。雖然目前RA仍無法治愈,但早期診斷和及時治療有助于減緩關節損害的進展,提高患者的生活質量[4]。RA的潛在發病機制尚未完全闡明,尋找與RA相關的發病關鍵基因并闡明其免疫機制有助于臨床提高疾病的早期診斷能力[5]。目前,基于機器學習算法的免疫浸潤和生物信息學分析可以挖掘新的診斷標志物,最近已應用于許多免疫相關疾病,包括潰瘍性結腸炎、白癜風、骨關節炎和銀屑病等[6]。鑒此,本研究通過生物信息學方法從RA數據集中篩選出差異表達基因(differentially expressed genes,DEGs),并進一步聯合機器學習算法尋找最佳的RA診斷基因,并對RA可能的發病機制及免疫浸潤機制進行探索,為RA的早期診斷以及靶向治療研發提供參考。

1 資料與方法

1.1數據資料檢索 在基因表達綜合(Gene Expression Omnibus,GEO) 數據庫(https://www.ncbi.nlm.nih.gov/geo/)以“rheumatoid arthritis”為檢索詞進行檢索。納入標準:(1)包含RA患者和健康對照數據;(2)樣本類別為關節滑膜組織;(3)芯片的平臺信息明確。排除標準:(1)芯片平臺上的基因探針ID無法轉化為基因名稱;(2)下載文件中的數據不完整或無效。根據上述納入、排除標準從GEO數據庫中選擇了3個基因數據集:GSE55235、GSE77298和GSE55457。將GSE55457作為獨立驗證數據集,將GSE55235和GSE77298作為聯合芯片訓練集。各芯片數據文件以及相應的平臺信息見表1。

表1 數據集信息

1.2數據預處理 使用“Limma”R包(版本)中的“normalizeBetweenArrays”函數來標準化GSE55235、GSE77298和GSE55457數據集,箱線圖用于規范化處理后的可視化。

1.3DEGs的篩選 使用R統計軟件中的Limma包研究聯合芯片訓練集的DEGs。篩選DEGs的截止值:adjustedP<0.05且|logFC|≥1。

1.4DEGs富集分析 應用R軟件的Cluster Profiler數據包對DEGs進行基因本體論(Gene Ontology,GO)富集分析和京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes,KEGG)富集分析,以P<0.05為標準篩選。

1.5機器學習篩選診斷基因 使用支持向量機-遞歸特征消除(support vector machines-recursive feature elimination,SVM-RFE)、最小絕對值收斂和選擇算子(least absolute shrinkage and selection operator,LASSO)和隨機森林(random forest,RF)三種機器學習算法篩選診斷基因,最后取三種方法獲得診斷基因的交集。

1.6診斷基因的驗證以及診斷效能分析 將GSE55457作為獨立的驗證數據集,采用成組t檢驗,以P<0.05判定差異有統計學意義,建立受試者工作特征(receiver operating characteristic,ROC)曲線,計算曲線下面積(area under the curve,AUC)值,評估研究指標診斷RA的效能。

1.7診斷模型的構建 整合診斷基因在驗證集數據集中的表達矩陣,使用logistic回歸分析構建診斷模型,應用R軟件構建可視化的列線圖來診斷RA患者,采用ROC曲線評估模型性能。

1.8免疫細胞浸潤分析 應用xCell算法[7]計算淋巴細胞、髓系細胞、基質細胞、干細胞以及其他免疫細胞的相對表達量,并采用秩和檢驗比較RA組與健康對照組(HC組)的表達差異,以P<0.05為差異有統計學意義。

2 結果

2.1數據預處理結果 應用R軟件Limma包對數據集GSE55235、GSE77298和GSE55457進行標化處理,批次校正后各數據集間的數據分布趨于一致。見圖1。

2.2DEGs分析結果 通過分析聯合芯片訓練集中的數據,總共篩選出RA組和HC組的DEGs共704個,其中上調476個,下調228個,并將結果可視化為火山圖(圖2?)和熱圖(圖2?)。

圖2 DEGs篩選結果圖

2.3DEGs的GO富集分析和KEGG富集分析結果 GO富集分析結果顯示,DEGs主要富集于白細胞介導的免疫、免疫應答的激活、白細胞遷移、淋巴細胞介導的免疫。見圖3?。KEGG富集分析結果顯示,DEGs主要富集于趨化因子信號通路、利什曼病、RA、金黃色葡萄球菌感染等。見圖3?。

圖3 DEGs的GO富集分析和KEGG富集分析結果圖

2.4機器學習篩選診斷基因結果 SVM-RFE算法挑選出27個診斷基因(見圖4??),構建LASSO回歸模型并進行交叉驗證,誤差最小值對應17個特征基因(見圖4??)。RF算法鑒定了39個特征基因(見圖4??)。取交集得到4個診斷基因:趨化因子CXC配體13(C-X-C motif chemokine ligand 13,CXCL13)、富含亮氨酸重復序列結構域15(leucine rich repeat containing 15,LRRC15)、多配體蛋白聚糖-1(syndecan 1,SDC-1)和核酸結合蛋白3(Y-box binding protein 3,YBX3)。見圖4?。

2.5診斷基因的驗證以及診斷效能分析 使用GSE55457芯片數據集對篩選出的4個診斷基因進行外部驗證,結果顯示,相對于HC組,YBX3基因表達在RA患者中顯著下調(P<0.05),而CXCL13、LRRC15和SDC-1基因表達在RA患者中顯著上調(P<0.05)。見圖5。這與訓練集中的基因表達趨勢相同。ROC曲線分析結果顯示,機器算法篩選得到的4個診斷基因在驗證數據集中對RA具有較高的診斷效能(AUC>0.8)。見圖6。

圖5 4個診斷基因在驗證數據集中表達水平比較圖

圖6 驗證數據集中4個診斷基因ROC曲線圖

2.6診斷RA的列線圖模型構建結果 基于驗證數據集的診斷基因表達矩陣,使用logistic回歸方法構建診斷預測模型,并通過列線圖對模型進行可視化。見圖7?。ROC曲線分析結果顯示,相對于任意一個單獨的診斷基因,組合模型對RA的診斷性能更高(AUC=0.985)。見圖7?。

圖7 診斷RA的列線圖模型構建結果圖

2.7免疫浸潤分析結果 與HC組相比,RA組B細胞、CD4+T細胞、樹突狀細胞和單核細胞水平顯著上調(P<0.05)。見圖8。

?淋巴細胞;?髓系細胞;?基質細胞;?干細胞和其他細胞。藍色為HC組,紅色為RA組;*P<0.05

3 討論

在本研究中,筆者對RA的基因表達芯片數據進行了生物信息學分析,共鑒定出了704個DEGs。GO和KEGG富集分析結果表明,DEGs在與免疫反應有關的生物過程中明顯富集。隨后聯合機器學習篩選出RA的診斷基因,最后得到4個診斷基因,即CXCL13、LRRC15、SDC-1和YBX3。在驗證數據集中,ROC曲線分析結果顯示這些基因對RA的診斷效能較高,具有較好的臨床應用前景。

3.1RA是一種以滑膜炎癥和進行性關節破壞為特征的慢性全身免疫介導疾病[8]。B細胞在RA發生中發揮了重要作用。自身反應性B細胞可產生直接參與關節損傷的IgG自身抗體[9]。CXCL13是B細胞一種關鍵的趨化因子,可通過調節局部免疫反應和抗體產生而參與疾病發生。研究表明,在RA患者中CXCL13水平顯著升高[10]。不僅如此,CXCL13還可反映RA疾病的嚴重程度。Meeuwisse等[11]研究表明,RA患者CXCL13基線水平越高,其關節破壞率越高。使用阿達木單抗控制RA病情后,CXCL13水平顯著下降,而在治療中斷時,CXCL13的水平上升至治療前水平[12]。動物實驗顯示,使用抗CXCL13的多克隆抗體治療可以降低小鼠關節疾病的嚴重程度并抑制生發中心的形成[13]。

3.2LRRC15目前被認為是一個有應用前景的抗癌靶點。而在RA方面,LRRC15在RA患者的滑膜組織中存在過表達,并與滑膜細胞的增殖、遷移、侵襲和血管生成能力增加以及促炎細胞因子的加速釋放有關。動物實驗顯示,關節炎的發展導致平衡型滑膜成纖維細胞萎縮,出現以LRRC15表達為標志的滑膜成纖維細胞輪廓,其功能是增強炎癥反應和基質分解過程[14]。敲低LRRC15表達不僅能夠在體外抑制滑膜細胞的侵襲性表型,而且在體內也能顯著抑制膠原誘導的關節炎小鼠的滑膜增殖,減少骨侵襲和破壞[2]。He等[5]也證實了LRRC15在RA模型大鼠滑膜組織中表達顯著增加。本研究結果也顯示LRRC15是參與RA微環境動態平衡失調的關鍵因子。

3.3SDC-1也被稱為CD138,它可與許多配體相互作用,引發與細胞黏附、血管生成、炎癥和組織修復相關的生物事件[15]。血清SDC-1被認為是炎癥活動的潛在標志物[16]。阻斷白介素(interleukin,IL)-34/SDC-1通路可減輕膠原誘導的關節炎癥和骨破壞,并加劇血管生成[17]。Deyab等[18]研究表明,RA患者血清C反應蛋白水平與SDC-1水平呈顯著正相關,當RA患者病情得到控制后,血清SDC-1水平顯著降低。在其他風濕性疾病中也觀察到相似的現象,如SDC-1水平與系統性紅斑狼瘡活動指數及抗dsDNA抗體水平呈正相關,活動性狼瘡腎炎患者的血清SDC-1水平也高于非活動性狼瘡腎炎患者及非腎炎患者[19]。

3.4YBX3是一個轉錄因子,參與調節上皮形態發生和穩態調節[20]。尤其在腫瘤性疾病方面的研究較為集中。YBX3通過調節膽汁酸生物合成途徑導致肝細胞癌轉移[21]。此外,YBX3高表達與直腸癌的深度浸潤相關,抑制YBX3表達可以減少體內腫瘤生長[22]。另外,抑制YBX3表達也可以提高結直腸癌細胞對化療藥物的敏感性[23]。然而,盡管YBX3的上調增加了癌細胞侵襲和腫瘤化療耐藥,但在某些腫瘤中也顯示出抗癌作用。有學者發現YBX3在低級別腎透明細胞癌中的表達水平較高,而在高級別腎透明細胞癌中表達降低[24]。本研究發現YBX3在RA中扮演重要角色,值得進一步通過基礎實驗開展研究。

3.5RA是一種全身炎癥性自身免疫性疾病,理想的治療策略應該是在組織損傷之前重新誘導自我耐受[25]。因此,明確疾病的免疫機制可為治療提供重要幫助。本研究通過xCell算法探討了RA的免疫特征,發現B細胞、CD4+T細胞、樹突狀細胞和單核細胞呈高表達,這可能是RA發生發展的關鍵免疫細胞。在RA中,B細胞主要向CD4+T輔助細胞呈遞自身抗原,外周血中的B細胞可分泌多種不同的細胞因子參與骨破壞,包括腫瘤壞死因子-α(tumor necrosis factor-α,TNF-α)、IL-6和IL-1β等[26]。目前,B細胞抑制劑利妥昔單抗在治療RA中取得良好效果?；そM織中的特殊成分和體內產生的內源性物質可由樹突狀細胞作為自身抗原呈現,激活CD4+T細胞并導致炎癥。Inamo等[27]認為可使用CD4+T細胞中的轉錄組數據對RA緩解和非緩解情況進行分類。T濾泡輔助細胞是CD4+T細胞的一種亞型,可以幫助B細胞調節抗體產生,從而進一步參與RA的發生[28]。在RA中,單核細胞離開血液并浸潤發炎的滑膜組織時,可以分化成外周血來源樹突狀細胞,通過促進IL-17的產生和向強效破骨細胞的轉化參與骨關節破壞;反之,用于控制RA炎癥的生物療法可調節單核細胞向樹突狀細胞的轉化過程。識別和控制外周血來源樹突狀細胞分化的環境介質以及潛在的分子信號通路,這可能是RA新療法開發的突破關鍵[29]。

綜上所述,本研究基于機器學習方法篩選出CXCL13、LRRC15、SDC-1和YBX3等對RA診斷具有應用前景的因子,闡釋了相關標志基因在RA中的生物學意義及可能的免疫機制,為RA的診斷和治療靶點開發提供了參考。但本研究由于數據資料的限制,無法評估生物標志物或免疫細胞與RA患者臨床特征的關聯性,還需要更多的臨床研究數據來進一步探索標志基因和相關免疫細胞在RA中的作用。