基于GEO數(shù)據(jù)庫(kù)分析高血壓腎病的差異表達(dá)基因及免疫細(xì)胞浸潤(rùn)

王益民 蘇小藝 蘇文革 彭偉 李焱 王怡斐

（1.山東中醫(yī)藥大學(xué)第一臨床醫(yī)學(xué)院，濟(jì)南 250014；2.山東中醫(yī)藥大學(xué)附屬醫(yī)院，濟(jì)南 250014）

高血壓腎病（hypertensive nephropathy，HN）由慢性高血壓引起，是終末期腎病的主要病因之一。高血壓是一種容易導(dǎo)致靶器官損傷的疾病，腎臟損傷是其重要的一部分，包括腎小球硬化和萎縮、足細(xì)胞丟失以及腎小管間質(zhì)纖維化等［1-2］。然而，有研究發(fā)現(xiàn)即使血壓降至目標(biāo)水平，也只能減緩而不能阻止HN的進(jìn)展［3］。腎小管間質(zhì)部分占腎臟總質(zhì)量的95%，間質(zhì)炎癥和纖維化是導(dǎo)致腎功能下降的重要原因，也是腎硬化癥的主要特征［4-5］。大量研究已證實(shí)免疫系統(tǒng)在高血壓的發(fā)病機(jī)制中發(fā)揮關(guān)鍵作用，淋巴細(xì)胞是高血壓發(fā)生和靶器官損傷的重要參與者［6］。目前尚無(wú)反映高血壓腎損傷的早期診斷標(biāo)志物能對(duì)HN的發(fā)生起到預(yù)警作用。基于此，本研究擬以腎小管間質(zhì)組織為研究對(duì)象，進(jìn)一步探討HN發(fā)病的相關(guān)機(jī)制，并發(fā)掘其潛在的診斷標(biāo)志物；并利用評(píng)估22種免疫細(xì)胞相對(duì)含量的CIBERSORT算法對(duì)HN腎小管間質(zhì)組織的免疫細(xì)胞浸潤(rùn)情況進(jìn)行分析研究，以揭示免疫細(xì)胞浸潤(rùn)模式。

1 資料與方法

1.1 GEO數(shù)據(jù)集的獲取 從GEO數(shù)據(jù)庫(kù)（https：//www.ncbi.nlm.nih.gov/geo/）下載HN患者和健康人（Normal）腎小管間質(zhì)組織樣本的數(shù)據(jù)集，GSE編號(hào)分別為GSE37455、GSE99325和GSE104954［4，7-8］。其中GSE37455和GSE104954均包括20例HN患者和21例健康人腎小管間質(zhì)組織樣本，GSE99325包括20例HN患者和4例健康人腎小管間質(zhì)組織樣本。所選數(shù)據(jù)集均具有完整的基因表達(dá)譜，同時(shí)不存在道德問(wèn)題、倫理問(wèn)題和其他利益沖突。本研究中實(shí)驗(yàn)組和驗(yàn)證組數(shù)據(jù)集選擇如表1所示。

表1 GEO數(shù)據(jù)集Tab.1 GEO datasets

1.2 數(shù)據(jù)處理和差異表達(dá)基因（differentially expressed genes，DEGs）的篩選 根據(jù)探針注釋文件將每個(gè)數(shù)據(jù)集中的探針名稱轉(zhuǎn)換為基因名稱，應(yīng)用R軟件（版本4.1.2）中“SVA”包的combat功能消除批次效應(yīng)，使用R軟件的“l(fā)imma”包對(duì)實(shí)驗(yàn)組中的40個(gè)HN樣本和24個(gè)Normal樣本進(jìn)行差異表達(dá)分析。調(diào)整后的假陽(yáng)性率P＜0.05和|log2FC|＞0.5的樣本作為DEGs的閾值。

1.3 功能富集分析 為進(jìn)一步探索DEGs的生物學(xué)意義，采用Metascape［9］在線分析平臺(tái)（https：//metascape.org/gp/index.html）對(duì)DEGs進(jìn)行功能富集分析。將P＜0.01、最小計(jì)數(shù)為3，且富集因子＞1.5的項(xiàng)分組到聚類中，相似性＞0.3的項(xiàng)由邊連接，并呈現(xiàn)為網(wǎng)絡(luò)圖。

1.4 機(jī)器學(xué)習(xí)篩選特征基因 本研究采用3種機(jī)器學(xué)習(xí)算法進(jìn)一步篩選DEGs。分別是最小絕對(duì)收縮和選擇算子（LASSO）、支持向量機(jī)遞歸特征消除法（SVM-RFE）和隨機(jī)森林（RF），在R軟件中分別用“glmnet”包、“e1071”包和“randomForest”包實(shí)現(xiàn)，然后將3種方法篩選的結(jié)果取交集，進(jìn)而確定HN的特征基因。

1.5 特征基因的診斷價(jià)值 分別使用實(shí)驗(yàn)組和驗(yàn)證組數(shù)據(jù)將得到的特征基因進(jìn)行驗(yàn)證，并繪制受試者工作特征（receiver operating characteristic，ROC）曲線，曲線下面積（area under curve，AUC）值用于確定特征基因在HN和Normal樣本中的診斷有效性。

1.6 免疫細(xì)胞浸潤(rùn)分析 將基因表達(dá)矩陣數(shù)據(jù)上傳至CIBERSORT，篩選出P＜0.05的樣本，得到免疫細(xì)胞浸潤(rùn)矩陣。然后利用R軟件中“corrplot”包繪制相關(guān)性熱圖，將22種免疫細(xì)胞浸潤(rùn)水平的相關(guān)性可視化；“ggplot2”包繪制小提琴圖，將免疫細(xì)胞浸潤(rùn)的差異可視化。

1.7 特征基因與免疫細(xì)胞浸潤(rùn)水平的相關(guān)性 為探討特征基因與免疫細(xì)胞浸潤(rùn)水平的相關(guān)性，采用R軟件進(jìn)行Spearman等級(jí)相關(guān)分析，并用“ggplot2”包進(jìn)行可視化。

2 結(jié)果

2.1 DEGs的結(jié)果 本研究對(duì)實(shí)驗(yàn)組40例HN和24例Normal樣本進(jìn)行了分析，在去除批次效應(yīng)后，使用“l(fā)imma”包對(duì)數(shù)據(jù)集進(jìn)行分析，共獲得277個(gè)DEGs，其中128個(gè)DEGs下調(diào)，149個(gè)DEGs上調(diào)，DEGs的結(jié)果顯示在熱圖和火山圖中（圖1、圖2）。

圖1 DEGs熱圖Fig.1 Heat map of DEGs

圖2 DEGs火山圖Fig.2 Volcano map of DEGs

2.2 功能富集分析結(jié)果 在Metascape在線分析平臺(tái)上進(jìn)行了DEGs的功能富集分析，圖3、圖4列出了排在前20的具有顯著意義的富集結(jié)果，并呈現(xiàn)出關(guān)系網(wǎng)絡(luò)，其中基因本體論（gene ontology，GO）生物過(guò)程主要與先天免疫反應(yīng)、對(duì)細(xì)胞因子的反應(yīng)、體液免疫反應(yīng)、細(xì)胞因子產(chǎn)生的正向調(diào)節(jié)、炎癥反應(yīng)和免疫系統(tǒng)過(guò)程的負(fù)向調(diào)節(jié)等有關(guān)。反應(yīng)組基因集（Reactome Gene Sets）主要與免疫系統(tǒng)中的細(xì)胞因子信號(hào)傳導(dǎo)、干擾素信號(hào)傳導(dǎo)和中性粒細(xì)胞脫顆粒等有關(guān)。維基通路（WikiPathways）主要與同種異體移植排斥等有關(guān)。

圖3 DEGs富集項(xiàng)的條形圖Fig.3 Bar graph of enriched terms of DEGs

圖4 DEGs富集項(xiàng)的網(wǎng)絡(luò)圖Fig.4 Network diagram of enrichment terms of DEGs

2.3 特征基因篩選結(jié)果 使用LASSO回歸算法確定17個(gè)基因作為HN的候選特征基因。使用SVMRFE算法確定40個(gè)基因作為HN的候選特征基因。使用RF算法得出最佳決策樹(shù)數(shù)目為136個(gè)，在DEGs中挑選出評(píng)分最高的30個(gè)基因作為候選特征基因。取三種不同算法得到基因的重疊部分。最后獲得了3個(gè)特征基因：CISH、GADD45A和ZFP36（圖5）。

圖5 三種機(jī)器學(xué)習(xí)算法篩選結(jié)果圖Fig.5 Screening results of three machine learning algorithms

2.4 特征基因的診斷價(jià)值結(jié)果 如圖6所示，3個(gè)特征基因在實(shí)驗(yàn)組和驗(yàn)證組中均表現(xiàn)出良好的診斷價(jià)值，在實(shí)驗(yàn)組數(shù)據(jù)集中，CISH的AUC值為0.963（95%CI：0.915～0.993），GADD45A的AUC值為0.953（95%CI：0.893～0.994），ZFP36的AUC值為0.916（95%CI：0.835～0.972）。在驗(yàn)證組數(shù)據(jù)集中，CISH的AUC值為0.793（95%CI：0.648～0.923），GADD45A的AUC值為0.775（95%CI：0.616～0.895），ZFP36的AUC值為0.755（95%CI：0.595～0.889）。

圖6 驗(yàn)證3個(gè)特征基因診斷價(jià)值的ROC曲線Fig.6 ROC curve for verifying diagnostic value of three characteristic genes

2.5 免疫細(xì)胞浸潤(rùn)分析結(jié)果 得到的22種免疫細(xì)胞浸潤(rùn)水平相關(guān)性熱圖如圖7，其中單核細(xì)胞和活化肥大細(xì)胞存在明顯的正相關(guān)關(guān)系（R=0.41），幼稚CD4+T細(xì)胞和活化肥大細(xì)胞存在明顯的正相關(guān)關(guān)系（R=0.37）。幼稚B細(xì)胞和記憶B細(xì)胞存在明顯的負(fù)相關(guān)關(guān)系（R=-0.64），活化肥大細(xì)胞和靜息肥大細(xì)胞存在明顯的負(fù)相關(guān)關(guān)系（R=-0.62），漿細(xì)胞和M2巨噬細(xì)胞存在明顯的負(fù)相關(guān)關(guān)系（R=-0.45）。免疫細(xì)胞浸潤(rùn)差異的小提琴圖顯示，與Normal組相比，HN組中腎小管間質(zhì)組織調(diào)節(jié)性T細(xì)胞和M1巨噬細(xì)胞浸潤(rùn)較多（圖8）。

圖7 HN中免疫細(xì)胞間的相關(guān)性分析熱圖Fig.7 Analysis of correlation between immune cells in HN by heat map

圖8 HN與Normal組間免疫細(xì)胞浸潤(rùn)差異的小提琴圖Fig.8 Violin diagram of the difference of immune cell infiltration between HN and Normal groups

2.6 特征基因與免疫細(xì)胞浸潤(rùn)水平的相關(guān)性 相關(guān)性分析結(jié)果顯示（圖9），CISH與靜息樹(shù)突狀細(xì)胞呈正相關(guān)（R=0.33，P=0.017），與M1巨噬細(xì)胞呈負(fù)相關(guān)（R=-0.41，P=0.002 6）。GADD45A與漿細(xì)胞（R=0.59，P=5.6e-06）、活化CD4+記憶T細(xì)胞（R=0.29，P=0.035）和幼稚CD4+T細(xì)胞（R=0.28，P=0.046）呈正相關(guān)，與M1巨噬細(xì)胞（R=-0.38，P=0.005）、M2巨噬細(xì)胞（R=-0.30，P=0.029）和調(diào)節(jié)性T細(xì)胞（R=-0.27，P=0.048）呈負(fù)相關(guān)。ZFP36與靜息樹(shù)突狀細(xì)胞（R=0.37，P=0.007 1）和活化肥大細(xì)胞（R=0.37，P=0.006 1）呈正相關(guān)，與M1巨噬細(xì)胞（R=-0.37，P=0.006 9）和靜息肥大細(xì)胞（R=-0.36，P=0.007 6）呈負(fù)相關(guān)。

圖9 3個(gè)特征基因與免疫細(xì)胞浸潤(rùn)水平的相關(guān)性Fig.9 Relationship between three characteristic genes and level of immune cell infiltration

3 討論

目前HN的治療仍以控制血壓為主，尚無(wú)特異性治療方式。由于人口老齡化和心血管疾病存活率的提高，預(yù)計(jì)HN的發(fā)病率在未來(lái)幾十年將進(jìn)一步增加，這也對(duì)HN的防治提出了新的挑戰(zhàn)。本研究通過(guò)分析包含40個(gè)HN和24個(gè)Normal腎小管間質(zhì)組織的基因表達(dá)譜，得到277個(gè)DEGs。利用Metascape在線分析平臺(tái)進(jìn)行生物信息學(xué)分析，發(fā)現(xiàn)DEGs與多種免疫相關(guān)過(guò)程有關(guān)，如先天免疫反應(yīng)、體液免疫反應(yīng)、對(duì)細(xì)胞因子的反應(yīng)、細(xì)胞因子產(chǎn)生的正向調(diào)節(jié)、免疫系統(tǒng)中的細(xì)胞因子信號(hào)傳導(dǎo)、干擾素信號(hào)傳導(dǎo)及中性粒細(xì)胞脫顆粒等有關(guān)。研究發(fā)現(xiàn)，先天免疫反應(yīng)的激活伴隨著炎癥反應(yīng)的發(fā)展，進(jìn)而引起纖維化、基質(zhì)沉積和進(jìn)行性腎損傷被認(rèn)為是包括高血壓腎損害在內(nèi)的多種腎臟疾病發(fā)病機(jī)制的關(guān)鍵因素［10］。體液免疫系統(tǒng)激活與高血壓發(fā)病之間的關(guān)聯(lián)也得到證實(shí)，系統(tǒng)性紅斑狼瘡（systemic lupus erythematosus，SLE）小鼠體液免疫系統(tǒng)激活和高血壓的發(fā)病具有明確的因果關(guān)系，用小鼠抗CD20抗體以消耗B細(xì)胞，可顯著減弱自身抗體的產(chǎn)生，并防止SLE小鼠模型高血壓的發(fā)展，強(qiáng)烈支持體液免疫系統(tǒng)激活可導(dǎo)致高血壓的發(fā)病［11-12］。多種細(xì)胞因子參與HN的發(fā)病機(jī)制，炎癥相關(guān)細(xì)胞因子如IL-1β、IL-6、單核細(xì)胞趨化蛋白-1（MCP-1）、細(xì)胞間黏附分子-1（ICAM-1）和TNF-α等參與了包括高血壓腎損傷在內(nèi)的炎癥反應(yīng)［13-14］。纖維化相關(guān)細(xì)胞因子如TGF-β、結(jié)締組織生長(zhǎng)因子（CTGF）等，在HN腎間質(zhì)纖維化的形成中發(fā)揮關(guān)鍵作用［15-17］。由此可見(jiàn)，免疫相關(guān)機(jī)制在HN的發(fā)病中扮演重要角色，為揭示HN的發(fā)病機(jī)制提供了新方向。

本研究采用三種機(jī)器學(xué)習(xí)算法進(jìn)一步篩選DEGs。LASSO是一種回歸分析算法，其使用正則化來(lái)提高預(yù)測(cè)精度［18］。支持向量機(jī)（SVM）是一種有監(jiān)督的機(jī)器學(xué)習(xí)技術(shù)，在分類和預(yù)測(cè)方面具有優(yōu)異的性能［19］。RF是通過(guò)集成學(xué)習(xí)的思想將多棵樹(shù)集成的一種算法，在處理多特征數(shù)據(jù)方面表現(xiàn)出優(yōu)異的性能［20］。這三種機(jī)器學(xué)習(xí)算法目前廣泛應(yīng)用于特征基因的篩選［21-22］。本研究基于這三種機(jī)器學(xué)習(xí)算法，最終在DEGs中篩選出3個(gè)特征基因，分別是CISH、GADD45A和ZFP36。細(xì)胞因子誘導(dǎo)的含SH2蛋白（CISH）屬于細(xì)胞因子信號(hào)傳導(dǎo)抑制因子（SOCS）家族成員，SOCS家族包括SOCS1～SOCS7和CISH，最新研究表明，SOCS家族成員可通過(guò)介導(dǎo)細(xì)胞因子信號(hào)傳導(dǎo)的負(fù)反饋抑制，在先天性和適應(yīng)性免疫應(yīng)答中起關(guān)鍵作用［23］。生長(zhǎng)停滯和DNA損傷誘導(dǎo)蛋白45α（GADD45A）是生長(zhǎng)停滯和DNA損傷誘導(dǎo)蛋白45（GADD45）蛋白家族成員，GADD45A在各種細(xì)胞功能中發(fā)揮至關(guān)重要的作用，包括DNA修復(fù)、細(xì)胞凋亡和DNA甲基化等［24］。鋅指蛋白36（ZFP36）是一種ARE結(jié)合蛋白，可促進(jìn)促炎細(xì)胞因子降解，例如TNF-α、粒細(xì)胞-巨噬細(xì)胞集落刺激因子（GM-CSF）、IL-6和環(huán)氧合酶-2（COX-2）等，在限制炎癥反應(yīng)方面具有重要作用［25］。3個(gè)特征基因在實(shí)驗(yàn)組和驗(yàn)證組均表現(xiàn)出良好的診斷價(jià)值，在生理功能上與細(xì)胞因子功能的調(diào)控、細(xì)胞功能的調(diào)節(jié)以及炎癥反應(yīng)的調(diào)節(jié)等相關(guān)，有望成為HN潛在的診斷標(biāo)志物和治療靶點(diǎn)。

免疫細(xì)胞浸潤(rùn)的結(jié)果顯示HN組中腎小管間質(zhì)組織調(diào)節(jié)性T細(xì)胞（regulatory T cells，Tregs）和M1巨噬細(xì)胞浸潤(rùn)較多，提示Tregs和M1巨噬細(xì)胞在HN中發(fā)揮關(guān)鍵作用。根據(jù)目前的研究，Tregs可影響免疫系統(tǒng)的多個(gè)環(huán)節(jié)來(lái)預(yù)防高血壓，以及減輕靶器官損傷［26］。缺乏Tregs會(huì)加重血管緊張素Ⅱ（Ang Ⅱ）依賴性高血壓，而反復(fù)過(guò)繼轉(zhuǎn)移Tregs則會(huì)減弱血壓升高［27］。Tregs產(chǎn)生的IL-10可通過(guò)減輕血管氧化應(yīng)激改善高血壓微血管內(nèi)皮功能［28］。Tregs還與補(bǔ)體系統(tǒng)相互作用。有研究表明，特異性靶向Tregs中的補(bǔ)體受體C3aR和C5aR可能是治療高血壓的另一種新方法［29］。M1巨噬細(xì)胞是極化巨噬細(xì)胞的一種，可產(chǎn)生活性氧和促炎細(xì)胞因子加劇炎癥，特別是TNF-α和IL-1β，介導(dǎo)了高血壓發(fā)病中的眾多機(jī)制［30］。還有研究指出，CD14+M1巨噬細(xì)胞通過(guò)其強(qiáng)烈表達(dá)血管緊張素轉(zhuǎn)換酶促進(jìn)血壓升高，提示它們可能通過(guò)RAS系統(tǒng)參與高血壓的發(fā)生［31］。22種免疫細(xì)胞浸潤(rùn)水平的相關(guān)性結(jié)果顯示，多種免疫細(xì)胞在HN的發(fā)病中具有相關(guān)性，單核細(xì)胞和活化的肥大細(xì)胞存在較強(qiáng)的正相關(guān)關(guān)系，幼稚B細(xì)胞和記憶B細(xì)胞存在較強(qiáng)的負(fù)相關(guān)關(guān)系。特征基因與免疫細(xì)胞浸潤(rùn)水平的相關(guān)性結(jié)果顯示，3個(gè)特征基因與M1巨噬細(xì)胞均具有顯著的負(fù)相關(guān)性，CISH和ZFP36與靜息樹(shù)突狀細(xì)胞均呈正相關(guān)，GADD45A與漿細(xì)胞、活化CD4+記憶T細(xì)胞等呈正相關(guān)。目前關(guān)于各免疫細(xì)胞之間、特征基因與免疫細(xì)胞之間相關(guān)性的基礎(chǔ)研究尚缺乏大規(guī)模驗(yàn)證，本研究結(jié)果具有一定的參考價(jià)值。

本研究也存在一定的局限性，數(shù)據(jù)來(lái)源依賴于GEO數(shù)據(jù)庫(kù)、樣本數(shù)量有限均可能造成分析結(jié)果的偏倚，同時(shí)還需要進(jìn)一步的實(shí)驗(yàn)驗(yàn)證。總之，本研究利用生物信息學(xué)和機(jī)器學(xué)習(xí)的方法對(duì)HN的DEGs進(jìn)行了分析研究，并探討了HN的免疫細(xì)胞浸潤(rùn)機(jī)制，對(duì)于HN發(fā)病機(jī)制的探討、診斷標(biāo)志物的篩選及治療均具有參考意義。