原發(fā)性干燥綜合征中鐵死亡相關生物標志物的生物信息學挖掘與診斷模型構建

燕紅梅, 吳舒婷, 巫燕琴, 董光富

(1. 華南理工大學 醫(yī)學院, 廣東 廣州, 510006;2. 南方醫(yī)科大學附屬廣東省人民醫(yī)院/廣東省醫(yī)學科學院 風濕免疫科, 廣東 廣州, 510080;3. 汕頭大學 醫(yī)學院, 廣東 汕頭, 515041; 4. 南方醫(yī)科大學, 廣東 廣州, 510515)

原發(fā)性干燥綜合征(pSS)是一種慢性自身免疫性疾病,主要臨床表現(xiàn)為口干和眼干,可伴發(fā)全身性損害[1]。目前, pSS的確切病因尚不明確,且癥狀特異性差,易誤診和延遲診斷,導致患者預后不佳。鐵死亡是一種新的程序性細胞死亡方式,在多種自身免疫性疾病的發(fā)病機制中起著關鍵作用。鐵死亡的特征是抗氧化劑谷胱甘肽消耗和谷胱甘肽過氧化物酶4水平降低,導致脂質過氧化物積累,從而引起廣泛的組織細胞死亡。異常死亡細胞釋放的組分具有潛在的免疫刺激效應,并可引發(fā)自身免疫反應[2]。研究[3-4]表明,鐵死亡參與了包括類風濕性關節(jié)炎(RA)和系統(tǒng)性紅斑狼瘡(SLE)在內的多種自身免疫性疾病的發(fā)生和發(fā)展。pSS作為自身免疫性疾病,其特點之一是持續(xù)的慢性炎癥,研究[5]表明, EB病毒感染等多種感染參與了pSS的發(fā)病和進展。鐵死亡是一種特殊類型的細胞死亡方式,可將慢性炎癥與氧化應激聯(lián)系起來[6]。當細胞被感染時,鐵死亡可以加劇炎癥反應,從而引起器官損害[7]。此外,氧化應激也是鐵死亡的主要特征之一,故鐵死亡可能在pSS的發(fā)病中扮演重要角色。本研究基于多種生物信息學手段識別與pSS相關的鐵死亡相關基因,通過機器學習算法去除冗余基因并構建診斷模型,分析鐵死亡相關基因在pSS發(fā)病過程中的可能作用,以期為pSS的診斷和治療提供新的視角。

1 資料與方法

1.1 數(shù)據(jù)來源

本研究從GEO數(shù)據(jù)庫(http://www.ncbi.nlm.nih.gov/geo/)中下載包含外周血單個核細胞(PBMC)樣本的5個RNAseq數(shù)據(jù)集進行分析,包括GSE51092(n=222)、GSE66795(n=164)、GSE84844(n=60)、GSE48378(n=27)和GSE208260(n=57)數(shù)據(jù)集,其中前2個為訓練集,后3個為驗證集。另從GEO數(shù)據(jù)庫下載1個單細胞數(shù)據(jù)集GSE157278, 包括5個來自pSS患者的樣本和5個來自健康供體的正常樣本。

1.2 鐵死亡相關基因收集

FerrDb是用于管理和鑒定鐵死亡相關標志物、調控因子以及鐵死亡關聯(lián)疾病的數(shù)據(jù)庫,本研究從FerrDb數(shù)據(jù)庫中獲取鐵死亡相關基因集,去除重復基因后,最終獲得484個基因。

1.3 數(shù)據(jù)處理和差異表達基因(DEGs)篩選

應用R軟件Combat包去除GSE51092數(shù)據(jù)集和GSE66795數(shù)據(jù)集的批次效應后進行合并,通過R軟件limma包進行差異分析,將調整后P值(P.adj)<0.05且|差異倍數(shù)(FC)|>1.3的基因作為DEGs。

1.4 基因本體論(GO)和京都基因與基因組百科全書(KEGG)富集分析

應用R軟件clusterProfiler包和org.Hs.eg.db包對pSS相關DEGs進行GO和KEGG富集分析,并應用ggplot2包將數(shù)據(jù)可視化。

1.5 加權基因共表達網(wǎng)絡分析(WGCNA)

基于R軟件WGCNA包構建pSS患者的加權基因共表達網(wǎng)絡,首先剔除離群的基因和樣本,確定軟閾值為4, 并構建無尺度網(wǎng)絡和拓撲重疊矩陣,然后將表達高度相關的基因劃分至同一基因模塊中,最后根據(jù)不同模塊與pSS的關聯(lián)性篩選出相關性最高的基因模塊。

1.6 診斷模型的構建和驗證

機器學習是人工智能的分支,本研究基于分類算法和回歸分析算法構建預測模型。最小絕對值收斂和選擇算子(LASSO)是一種正則化模型,通常與線性分類器算法支持向量機(SVM)同時使用。LASSO通過對高維數(shù)據(jù)集中無關變量的系數(shù)進行懲罰或縮減,最終獲得1個變量較少的模型,具有從高度相關的變量中選擇相關特征的能力[8]。SVM是一種監(jiān)督學習算法,通過建立2個類別之間的決策邊界,可從1個或多個特征向量中獲得預測標簽,結合多個機器學習方法能進一步提高模型的預測能力。

從DEGs、鐵死亡相關基因和WGCNA紅色模塊中篩選出關鍵基因,通過LASSO和SVM-遞歸特征消除(RFE)這2種機器學習算法去除冗余基因。繪制受試者工作特征(ROC)曲線,評估關鍵基因的診斷效能,將曲線下面積(AUC)大于0.8的關鍵基因重新輸入LASSO模型中構建最終的診斷模型。將GSE84844、GSE48378、GSE208260數(shù)據(jù)集作為驗證隊列,基于ROC曲線評估診斷模型的診斷性能。

1.7 免疫浸潤分析

采用CIBERSORT算法,估計22種不同免疫細胞亞型的浸潤情況。從CIBERSORT網(wǎng)站(http://cibersort.stanford.edu/)獲取LM22文件,通過R軟件對結果進行可視化處理,應用ggcorrplot包可視化表示關鍵基因與免疫細胞的相關性。

1.8 單細胞RNA測序的質控、降維和注釋

使用R軟件Seurat包進行單細胞測序數(shù)據(jù)分析,過濾表達少于3個基因或少于200個基因的細胞,并將表達超過2 500個基因或線粒體基因超過10%的細胞排除。通過主成分分析(PCA)對單細胞測序數(shù)據(jù)進行降維處理,應用harmony軟件包合并單細胞樣本,通過FindNeighbors和FindClusters(Dim=30, Resolution=1.9)進行細胞聚類,然后用統(tǒng)一流形逼近和投影(UMAP)圖展示聚類結果。應用FindAllMarkers函數(shù)確定每個聚類的標記基因。將主成分進行聚類和UMAP表示,最終確定聚類。基于既往研究報道的細胞標記基因,對每個聚類進行手動注釋,最終注釋出細胞類型。

1.9 不同PBMC細胞類型中關鍵基因的表達情況

通過ggplot2和ggalluvial包繪制每個樣本中不同細胞類型的百分比,通過plot1cell包展示PBMC中各種細胞類型關鍵基因的表達情況。

1.10 統(tǒng)計學分析

應用R軟件(版本4.2.1)進行統(tǒng)計學分析,采用Wilcoxon秩和檢驗或學生t檢驗分析2組間的差異,采用Spearman秩相關檢驗確定變量之間的相關性。所有統(tǒng)計檢驗采用雙側檢驗,P<0.05為差異有統(tǒng)計學意義。

2 結 果

2.1 pSS中DEGs的鑒定及其GO和KEGG富集分析

本研究將2個大樣本數(shù)據(jù)集(GSE51092、GSE66795數(shù)據(jù)集)作為訓練集,應用R軟件Combat包消除兩者的批次效應后進行合并。去除批次效應后, UMAP圖顯示2個數(shù)據(jù)集的樣本交織在一起,提示批次效應已成功消除,見圖1A、圖1B。合并后的數(shù)據(jù)集中共鑒定出265個DEGs(P.adj<0.05且|FC|>1.3), 其中179個上調, 86個下調。

應用R軟件對265個DEGs進行GO和KEGG富集分析,結果見圖1C、圖1D。GO富集分析包括生物學過程(BP)、細胞成分(CC)和分子功能(MF), DEGs的BP包括防御病毒反應、防御共生反應、對病毒的反應等, CC主要富集于ISGF3復合物、含有膠原的細胞外基質和TAP復合物等, MF主要包括雙鏈R結合、D+蛋白質ADP核糖轉移酶活性、MHC I類b蛋白質結合和D+ ADP核糖轉移酶活性; KEGG富集分析顯示, DEGs參與多種信號通路,包括流感A病毒、丙型肝炎病毒、Epstein-Barr(EB)病毒感染等。

A: 去除批次效應前UMAP圖; B: 去除批次效應后UMAP圖; C: DEGs的GO富集分析結果; D: DEGs的KEGG通路富集分析結果。圖1 pSS相關DEGs的功能分析結果

2.2 WGCNA鑒定與pSS相關的模塊

當軟閾值確定為β=4,R2=0.86, 網(wǎng)絡符合無標度網(wǎng)絡的分布,見圖2A、圖2B。應用混合動態(tài)剪切樹法識別相似的模塊并進行合并,最終獲得27個基因模塊,見圖2C。模塊-臨床特征相關性分析結果顯示,紅色模塊與pSS的相關性最高(r=0.44,P<0.01), 見圖2D。紅色模塊內基因成員與疾病狀態(tài)高度相關(r=0.82), 見圖2E。將紅色模塊基因、鐵死亡相關基因、DEGs取交集,共得到8個關鍵基因,見圖2F。

A、B: 基于比例獨立性、平均連接分析選擇β=4作為軟閾值構建網(wǎng)絡; C: 基因樹顯示不同基因共表達模塊(各模塊以不同的唯一顏色表示);D: 熱圖顯示基因模塊與pSS的相關性(紅色模塊與pSS顯著相關); E: 紅色模塊內基因成員身份與基因重要性的相關性分析散點圖; F: 紅色模塊基因、鐵死亡相關基因、DEGs取交集的韋恩圖。圖2 基于WGCNA鑒定與pSS相關的模塊

2.3 與鐵死亡相關的關鍵基因的鑒定及pSS診斷模型建立

為了去除冗余基因,應用LASSO模型篩選出4個與鐵死亡相關的診斷pSS的關鍵基因,見圖3A; 為了簡化診斷模型并提高診斷模型的準確性,基于SVM-RFE算法對取交集得到的8個關鍵基因進行二次篩選,見圖3B。2種機器學習方法篩選出4個共同的基因,即TRIM21、PARP9、PARP12和PARP14。繪制ROC曲線評估這4個基因的診斷價值(圖3C～圖3F), 將AUC>0.8的關鍵基因PARP9、PARP12、PARP14作為pSS的生物標志物。將3個關鍵基因輸入LASSO模型中,得到最終模型方程,即風險評分=-0.523 290 132 2+0.000 010 735 4×PARP12+0.115 903 498 1×PARP9+0.027 443 175 8×PARP14。

A: 基于LASSO模型篩選pSS診斷生物標志物(最佳基因數(shù)為4個,見曲線最低點); B: 基于SVM-RFE算法的關鍵基因圖; C、D、E、F: 訓練集中PARP9、PARP12、PARP14、TRIM21診斷pSS的ROC曲線。圖3 基于機器學習算法鑒定與鐵死亡相關的關鍵基因

2.4 診斷模型的驗證結果

為了驗證診斷模型的準確性,本研究將GSE84844、GSE208260、GSE48378數(shù)據(jù)集作為驗證集。3個驗證集中, 3個關鍵基因的表達情況見圖4A～圖4C; 3個關鍵基因在pSS患者(pSS組)和健康個體(對照組)中的表達水平存在差異,見圖4D～圖4F。ROC曲線顯示,診斷模型在GSE84844、GSE208260、GSE48378數(shù)據(jù)集中的AUC分別為0.848、0.853、1.000, 見圖4G～圖4I, 表明診斷模型具有區(qū)分pSS患者和健康個體的能力。

A: 驗證集GSE84844中3個關鍵基因的表達情況; B: 驗證集GSE208260中3個關鍵基因的表達情況; C: 驗證集GSE48378中3個關鍵基因的表達情況; D、E、F: 在GSE84844、GSE208260、GSE48378數(shù)據(jù)集中進行的診斷模型分類; G、H、I: 診斷模型在GSE84844、GSE208260 GSE48378數(shù)據(jù)集中的ROC曲線。圖A、B、C中,兩者比較, ??P<0.01, ???P<0.001, ????P<0.000 1。圖4 診斷模型的驗證結果

2.5 免疫細胞浸潤分析

本研究采用直方圖顯示每個樣本免疫細胞分布的總體情況,結果見圖5A。箱線圖表明,相較于對照組, pSS組患者的活化樹突狀細胞(DC)、幼稚B細胞、活化的CD4+記憶型T細胞水平較高,而靜止期CD4+記憶型T細胞和記憶型B細胞水平較低,見圖5B。22種免疫細胞之間的相關性分析表明, CD8+T細胞與中性粒細胞呈負相關(r=-0.46), 調節(jié)性T細胞(Tregs)與記憶型B細胞呈正相關(r=0.39), 見圖5C。進一步分析3個關鍵基因與22種免疫細胞的相關性,結果顯示, 3個關鍵基因(PARP12、PARP14、PARP9)均與活化DC呈正相關(r=0.76、0.75、0.72), 與M0型巨噬細胞呈負相關(r=-0.24、-0.22、-0.17), 見圖5D。

A: 每個樣本中22種免疫細胞的分布直方圖; B: 健康人和pSS患者之間22種免疫細胞的相關性分析結果(兩者比較, ??P<0.01, ????P<0.000 1, ns P>0.05); C: pSS組和非pSS組的免疫浸潤箱線圖; D: 3個關鍵基因與免疫細胞及免疫功能的相關性分析結果。圖5 免疫細胞浸潤分析

2.6 3個關鍵基因在pSS患者PBMC中的表達

為了進一步探討3個關鍵基因在pSS患者PBMC中的表達情況,本研究基于pSS單細胞數(shù)據(jù)集進一步分析。該數(shù)據(jù)集共有10個樣本,即5個pSS患者的樣本(pSS組)和5個健康供體的對照樣本(正常組),包含21 176個基因和54 152個細胞,其中pSS組包含28 938個細胞,正常組包含25 214個細胞。應用Seurat包中的UMAP方法將這些細胞分為31個聚類,根據(jù)既往文獻報道的細胞標志物,手動注釋出9種細胞類型,見圖6A、圖6B。正常組和pSS組各細胞類型中3個關鍵基因的表達情況見圖6C～圖6E, 結果顯示,PARP9和PARP14在pSS組的巨噬細胞、常規(guī)DC(cDC)、CD14+單核細胞中高度表達,PARP12在巨噬細胞、cDC、CD14+單核細胞中的表達水平未發(fā)生變化,但pSS組表達這3個基因的細胞比例增加。

A: 點圖顯示每個聚類的細胞標記基因的平均表達水平; B: 來自5個pSS患者和5個健康供體的單細胞UMAP圖; C、D、E: 各種細胞類型中3個中心基因的表達。圖6 pSS單細胞RNA測序譜

3 討 論

pSS的癥狀及診斷指標特異性差,常與其他自身免疫性疾病(如SLE和RA)重疊[9], 易導致診斷延遲和誤診,故尋找新的特異性生物標志物有助于pSS的診斷和治療。本研究共鑒定出265個DEGs, GO分析結果表明其與響應病毒感染密切相關。既往研究[10]表明,病毒可通過各種機制引起自身免疫性疾病,包括旁觀者活化、分子模擬、表位擴散和B細胞免疫活化等。KEGG通路分析同樣表明, DEGs富集于與病毒相關的通路,包括流感病毒A和EB病毒感染。流感病毒A可通過其NS1mRNA的另類閱讀框架觸發(fā)強烈的CD8+T細胞反應或刺激漿細胞樣樹突狀細胞產(chǎn)生α干擾素(IFN-α)而誘導自身免疫性。另有研究[11]表明, EB病毒也參與了pSS的發(fā)病機制。因此,研究人員或有必要更深入地探討病毒引起的自身免疫性在pSS中的作用。

本研究運用多種生物信息學方法識別pSS中與鐵死亡相關的關鍵基因,并分析了3個與鐵死亡相關的關鍵基因(PARP9、PARP12和PARP14)對pSS發(fā)病的診斷價值。鐵死亡是一種新的細胞死亡方式[12], 研究[4, 13-14]發(fā)現(xiàn),鐵死亡在多種自身免疫性疾病(如SLE、RA和癌癥)的發(fā)生和進展中具有重要作用。既往研究[15]表明,干擾PARP表達可以調節(jié)溶質載體家族7成員11(SLC7A11), 并參與乳腺癌易感基因(BRCA)介導的卵巢癌中的鐵死亡過程。在RA中,PARP9是一種具有不同甲基化位點的基因,可影響T淋巴細胞的轉錄表達,表明PARP9在RA的發(fā)病機制中具有重要性[16]。JIANG Z H等[17]發(fā)現(xiàn),在SLE患者的多種免疫細胞中,PARP12表現(xiàn)出低甲基化水平。另有研究[18]表明,基于PARP14、ATP10A和MX1這3個基因構建的評分模型在抗Ro 60陽性的自身免疫性疾病患者中的評分始終高于抗Ro 60陰性的患者。以上研究表明,PARP9、PARP12和PARP14在自身免疫性疾病中可能發(fā)揮重要作用,然而PARP9、PARP12和PARP14在pSS發(fā)展中的作用尚不明確。

本研究免疫浸潤分析結果表明, pSS患者存在活化DC上調, B細胞與T細胞亞群比例失調。進一步觀察pSS患者PBMCs的單細胞數(shù)據(jù)集,本研究發(fā)現(xiàn)DC、B細胞、單核細胞和巨噬細胞中PARP9和PARP14上調,表達PARP12的細胞在B細胞和巨噬細胞中的比例略微增加。DC在維持免疫穩(wěn)態(tài)和自身耐受性方面發(fā)揮著重要作用,研究[19]已證實cDC破壞CD4+T細胞自身耐受性,導致自身抗體產(chǎn)生、脾腫大和Th1、Th17效應細胞的增加。一項研究[20]表明,在pSS患者中, IFN-α可高度誘導DC中PARP9的表達,且DC的抗原攝取、處理和呈遞功能均發(fā)生改變,這增加了Ⅰ型IFN的產(chǎn)生并誘導CD4+T細胞擴增,促進了pSS的發(fā)展。此外,在RA中,PARP9的甲基化與T細胞轉錄表達的增加相關[16]。在白細胞介素(IL)-4存在的情況下,PARP14介導信號轉導與轉錄激活因子6(Stat6)與其靶基因結合,促進B細胞失調,這是pSS的病理標志之一[21]。同時,PARP14通過與Stat6相互作用調控細胞因子IL-4、IL-5、IL-13的表達,從而促進Th9細胞的發(fā)生發(fā)展和Th2細胞的分化[22]。這些研究表明,PARP家族可能通過異常調控免疫細胞參與pSS的發(fā)病。

本研究通過差異基因分析和WGCNA鑒定出與pSS相關的鐵死亡基因,并應用機器學習算法去除冗余基因,篩選出3個可靠的生物標志物(PARP9、PARP12和PARP14)用于構建診斷模型; 在3個驗證集中,該診斷模型均表現(xiàn)出高診斷效能,其AUC分別為0.848、0.853和1.000; 免疫浸潤分析發(fā)現(xiàn), pSS患者活化DC水平較高, T細胞與B細胞比例失衡, 3個關鍵基因與活化DC均呈正相關; 基于單細胞數(shù)據(jù)集驗證了3個關鍵基因在免疫細胞中的表達情況。綜上所述,本研究基于生物信息學手段識別出pSS中與鐵死亡相關的3個關鍵基因并構建pSS診斷模型,或可為pSS的診斷提供新工具。