基于加權(quán)基因共表達(dá)網(wǎng)絡(luò)鑒定腎移植術(shù)后排斥反應(yīng)中巨噬細(xì)胞M1亞型相關(guān)基因

董博清 李楊 石玉婷 張靜 馮新順 鄭瑾 李瀟 丁小明 薛武軍

與透析相比較，腎移植能明顯改善終末期腎病患者的生存質(zhì)量，并且減少遠(yuǎn)期并發(fā)癥[1]。隨著腎移植術(shù)后新型免疫抑制藥的應(yīng)用，1年移植腎存活率達(dá)95%，然而移植腎長期存活率依舊沒有明顯改善[2]。因此需要明確影響移植物長期存活的危險因素，對具有高危危險因素受者進(jìn)行早期干預(yù)以改善遠(yuǎn)期結(jié)局。

目前，排斥反應(yīng)仍然是影響移植腎存活的最主要原因[3]。在排斥反應(yīng)期間，多種細(xì)胞浸潤移植腎，目前的研究主要集中在適應(yīng)性免疫細(xì)胞，然而，越來越多的證據(jù)表明固有免疫中的巨噬細(xì)胞與移植腎功能和臨床預(yù)后密切相關(guān)[4-5]。巨噬細(xì)胞屬于單核吞噬細(xì)胞系統(tǒng)，是具有高度可塑性的細(xì)胞[6]。經(jīng)典的巨噬細(xì)胞M1亞型，由干擾素（interferon，IFN）-γ，脂多糖和腫瘤壞死因子（tumor necrosis factor，TNF）-α誘導(dǎo)，屬于巨噬細(xì)胞亞群中促炎細(xì)胞亞群，已被證明與移植物預(yù)后不良密切相關(guān)[7]，然而由于巨噬細(xì)胞極化的時間以及空間特性，目前對于調(diào)控M1亞型極化及募集在移植腎中的具體機(jī)制尚不清楚。因此，有必要鑒定排斥反應(yīng)中參與調(diào)控巨噬細(xì)胞M1亞型的基因，更加深入了解這些基因如何參與移植腎損傷并影響預(yù)后，開發(fā)診斷及預(yù)后的生物標(biāo)志物。本研究通過加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA）以及差異表達(dá)分析確定了腎移植術(shù)后排斥反應(yīng)中M1亞型相關(guān)的基因，并基于這些基因構(gòu)建了預(yù)測移植腎存活的風(fēng)險模型。

1 資料與方法

1.1 數(shù)據(jù)集處理

從基因表達(dá)綜合（Gene Expression Omnibus，GEO）數(shù)據(jù)庫下載腎移植術(shù)后的微陣列數(shù)據(jù)集GSE36059以及GSE21374[8-9]。兩個數(shù)據(jù)集均基于GPL570平臺，GSE36059數(shù)據(jù)集中包含排斥反應(yīng)和穩(wěn)定移植物樣本，使用該數(shù)據(jù)集篩選巨噬細(xì)胞M1亞型相關(guān)的基因；GSE21374數(shù)據(jù)集中包含了移植物丟失的隨訪數(shù)據(jù)，使用該數(shù)據(jù)集構(gòu)建風(fēng)險模型并進(jìn)行驗證。使用R語言limma包校正兩個數(shù)據(jù)集[10]。

1.2 免疫浸潤及加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

CIBERSORT是一種使用基因表達(dá)數(shù)據(jù)估計混合細(xì)胞群中成員細(xì)胞類型的豐度的反卷積算法。對GSE36059數(shù)據(jù)集進(jìn)行CIBERSORT反卷積分析，可計算出各個樣本中免疫細(xì)胞的含量。WGCNA是用來描述基因關(guān)聯(lián)模式的系統(tǒng)生物學(xué)方法，可以用來鑒定高度協(xié)同變化的基因集，并根據(jù)基因集的內(nèi)連性和基因集與表型之間的關(guān)聯(lián)鑒定候補(bǔ)生物標(biāo)記基因。為了篩選出具有共表達(dá)性質(zhì)的基因，使用GSE36059數(shù)據(jù)集中均絕對誤差（mean absolute deviation，MAD）前5 000的基因構(gòu)建共表達(dá)網(wǎng)絡(luò)。首先進(jìn)行數(shù)據(jù)質(zhì)控，剔除離群樣本。然后根據(jù)無標(biāo)度拓?fù)鋽M合指數(shù)R2以及平均連通性確定最佳的軟閾值。通過動態(tài)切割樹算法識別模塊，設(shè)置模塊內(nèi)最小基因數(shù)為50，將構(gòu)建的模塊聚類以合并相似模塊。最后進(jìn)行相關(guān)性分析，分析模塊和免疫細(xì)胞的相關(guān)性，選擇和M1亞型相關(guān)系數(shù)r最高的模塊進(jìn)行后續(xù)分析。在模塊中選擇和模塊的|r|>0.8且和M1亞型的|r|>0.4的基因作為M1亞型相關(guān)基因。

1.3 差異表達(dá)分析

在GSE36059數(shù)據(jù)集中使用limma包進(jìn)行差異表達(dá)分析，篩選出排斥反應(yīng)樣本和穩(wěn)定移植物中差異表達(dá)基因（differentially expressed gene，DEG）。差異倍數(shù)|（Log2 fold change，Log2FC）|>1且校正后P<0.05的基因為DEG。

1.4 基于差異表達(dá)的M1亞型相關(guān)基因的預(yù)后模型構(gòu)建

將WGCNA以及差異表達(dá)分析得到的基因取交集，最終篩選出排斥反應(yīng)樣本和對照樣本中差異表達(dá)的M1-DEG。隨后，按照7∶3拆分GSE21374數(shù)據(jù)集為訓(xùn)練隊列以及驗證隊列。在訓(xùn)練隊列中，使用M1-DEG進(jìn)行單因素Cox分析，篩選P<0.05的基因構(gòu)建模型。最小絕對收縮和選擇算法（least absolute shrinkage and selection operator，LASSO），是通過生成一個懲罰函數(shù)以壓縮回歸模型中變量的系數(shù)，解決嚴(yán)重共線性的問題[11]。將上一步篩選出的基因帶入LASSOCox回歸中，根據(jù)最小均方誤差（minimum mean square error，MMSE）選擇懲罰系數(shù)λ進(jìn)一步篩選構(gòu)建風(fēng)險模型的變量，最后生成列線圖可視化模型。

1.5 模型評估

隨后通過受試者工作特征曲線下面積（area under curve，AUC）評估基于M1-DEG的風(fēng)險模型預(yù)測1年及3年移植物存活的能力。根據(jù)風(fēng)險評分的中位數(shù)將患者分為高風(fēng)險組和低風(fēng)險組，使用對數(shù)秩檢驗（log-rank）明確不同組移植物存活是否存在差異。

1.6 免疫浸潤及人類白細(xì)胞抗原相關(guān)基因分析

使用CIBERSORT分析高風(fēng)險組和低風(fēng)險組22種浸潤的免疫細(xì)胞，使用Wilcoxn秩和檢驗分析組間差異，P<0.05為差異有統(tǒng)計學(xué)意義。使用熱圖描述高風(fēng)險和低風(fēng)險組人類白細(xì)胞抗原（human leukocyte antigen，HLA）相關(guān)基因的分布。

1.7 基因集富集分析

通過基因集富集分析（gene set enrichment analysis，GSEA），進(jìn)一步闡明高風(fēng)險組富集的生物學(xué)過程以及京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路。使用enrichplot包可視化GSEA結(jié)果展示標(biāo)準(zhǔn)化富集分?jǐn)?shù)（normalized enrichment score，NES）前5的生物學(xué)過程及KEGG通路，P<0.05為差異有統(tǒng)計學(xué)意義。

1.8 預(yù)測調(diào)控預(yù)后基因的miRNA

“multiMiR”包集成了7個在線數(shù)據(jù)庫的微小核糖核酸（microRNA，miRNA，miR）和預(yù)后基因之間相互作用的信息[12]。為了確保篩選出可靠的和預(yù)后基因互作的miRNA，僅納入出現(xiàn)在≥2個數(shù)據(jù)庫且通過實驗證實的相互作用關(guān)系。最后，使用可視化軟件Cytoscape繪制miRNA-預(yù)后基因間的調(diào)控關(guān)系。

1.9 統(tǒng)計學(xué)方法

所有統(tǒng)計分析基于R語言（4.1.2版本）?；蚰K與臨床特征相關(guān)性采用Pearson相關(guān)性分析。采用單因素Cox分析篩選與移植物存活密切相關(guān)的基因，并采用LASSO-Cox回歸分析選擇懲罰系數(shù)λ進(jìn)一步篩選構(gòu)建風(fēng)險模型的變量。采用log-rank檢驗分析移植物存活率。采用受試者工作特征曲線評估1、3年移植物存活。雙側(cè)P<0.05為差異有統(tǒng)計學(xué)意義。

2 結(jié) 果

2.1 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

首先進(jìn)行數(shù)據(jù)質(zhì)控，剔除GSE36059數(shù)據(jù)集中離群樣本（圖1A）。根據(jù)R2以及平均連通性選擇5作為軟閾值構(gòu)建模塊（圖1B）。通過動態(tài)剪切算法最初識別出10個模塊，而后合并相似的棕色以及黃色模塊（圖1C），最終得到9個模塊（圖1D）。將CIBERSORT計算出的GSE36059各個樣本中浸潤的免疫細(xì)胞的含量和模塊進(jìn)行Pearson相關(guān)性分析，其中棕色模塊和巨噬細(xì)胞M1亞型相關(guān)性最高（圖1E、F，r=0.45，P<0.001）。最后，篩選出棕色模塊中54個和模塊及M1亞型相關(guān)性高的基因作為M1亞型相關(guān)基因（圖1G）。

圖1 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析Figure 1 Analysis of weighted gene co-expression networks

2.2 確定差異表達(dá)的M1亞型相關(guān)基因

根據(jù)在GSE36059數(shù)據(jù)集中篩選出42個上調(diào)DEG（圖2A）。將DEG和WGCNA得到的M1亞型相關(guān)基因取交集后，最終得到14個M1-DEG（圖2B）。

圖2 差異表達(dá)的火山圖和韋恩圖Figure 2 Differentially expressed volcano plots and Wayne diagrams

2.3 基于M1-DEG的風(fēng)險模型構(gòu)建

將篩選出的14個M1-DEG帶入GSE21374訓(xùn)練集中，使用單因素Cox分析進(jìn)一步篩選與移植物存活密切相關(guān)的基因，14個M1-DEGP<0.05。而后使用LASSO-Cox回歸（圖3A），根據(jù)最小均方誤差選擇λ（圖3B），最終篩選出Toll樣受體8（Toll-like receptor 8，TLR8）、Fc γ 受 體 1B（Fc gamma receptor 1B，F(xiàn)CGR1B）、BCL2相關(guān)蛋白A1（BCL2 related protein A1，BCL2A1）、組織蛋白酶S（cathepsin S，CTSS）、鳥苷酸結(jié)合蛋白2（guanylate binding protein 2，GBP2）及半胱氨酸天冬氨酸蛋白酶招募域家族成員16（caspase recruitment domain family member 16，CARD16）構(gòu)建風(fēng)險預(yù)后模型[風(fēng)險 評 分 =?0.921 9×TLR8+（?0.591 7）×FCGR1B+（?1.352 2）×BCL2A1+1.401 7×CTSS+2.056 6×GBP2+0.904 2×CARD16]。使用列線圖可視化風(fēng)險模型預(yù)測1年及3年移植物存活，風(fēng)險評分= ?0.921 9×TLR8+（?0.591 7）×FCGR1B+（?1.352 2）×BCL2A1+1.401 7×CTSS+2.056 6×GBP2+0.904 2×CARD16（圖3C）。

圖3 基于M1-DEG構(gòu)建風(fēng)險預(yù)后模型Figure 3 Risk prognosis model based on M1-DEG

2.4 模型評估

根據(jù)公式計算各個受者的風(fēng)險評分，根據(jù)風(fēng)險評分中位數(shù)將受者分為高風(fēng)險組和低風(fēng)險組（圖4A、圖5A）。在訓(xùn)練集中，風(fēng)險評分高的受者移植物存活時間短（圖4B），基于6個M1-DEG的風(fēng)險模型預(yù)測1年及3年移植物存活的AUC分別為0.918和0.877（圖4C）。生存分析表明，高風(fēng)險組移植物存活率低于低風(fēng)險組（圖4D，P<0.001）。在驗證集中的結(jié)果和訓(xùn)練集相近（圖5B、D），模型預(yù)測移植物1年及3年存活的AUC分別為0.765及0.736（圖5C）。

圖4 6個基因模型在訓(xùn)練集的表現(xiàn)Figure 4 Six gene models performed in the training set

圖5 6個基因模型在驗證集的表現(xiàn)Figure 5 Six gene models performed in the validation set

2.5 免疫浸潤分析及HLA相關(guān)基因分析

在高風(fēng)險組中，靜息及活化的CD4+記憶T細(xì)胞、γδT細(xì)胞、巨噬細(xì)胞M1亞型浸潤；而記憶B細(xì)胞、巨噬細(xì)胞M0亞型、巨噬細(xì)胞M2亞型減少（均為P<0.05，圖6A）。在高風(fēng)險組中，HLAⅠ類基因HLA-A、B、C表達(dá)上調(diào)，部分HLAⅡ類位點(diǎn)的表達(dá)上調(diào)（圖6B）。

圖6 高風(fēng)險組和低風(fēng)險組的免疫浸潤及相關(guān)基因分布分析Figure 6 Immune infiltration and related gene distribution analysis between high-risk group and low-risk group

2.6 基因集富集分析

GSEA表明高風(fēng)險組中明顯富集的生物學(xué)過程包括免疫反應(yīng)的活化、適應(yīng)性免疫反應(yīng)、基于免疫球蛋白超家族結(jié)構(gòu)域構(gòu)建的免疫受體體細(xì)胞重組的適應(yīng)性免疫反應(yīng)、α/β T細(xì)胞的活化、α/β T細(xì)胞的分化（圖7A）；高風(fēng)險組富集的KEGG通路包括移植物排斥反應(yīng)、抗原提呈和加工、自身免疫性甲狀腺疾病、細(xì)胞黏附分子、趨化因子信號通路（圖7B）。

圖7 基因集富集分析Figure 7 Analysis of gene set enrichment

2.7 調(diào)控預(yù)后基因的miRNA預(yù)測

基于7個在線數(shù)據(jù)庫，最后確認(rèn)了13個miRNA和4個預(yù)后基因間的14個相互作用關(guān)系，其中CTSS與8個miRNA相互作用、BCL2A1和GBP2與3個miRNA相互作用、FCGR1B與1個miRNA相互作用（圖8）。

圖8 miRNA和預(yù)后基因互作網(wǎng)絡(luò)Figure 8 miRNA and prognostic gene interaction network

3 討 論

盡管免疫誘導(dǎo)以及以鈣調(diào)磷酸酶抑制劑為核心的三聯(lián)免疫抑制方案提高了移植腎的早期存活率，但移植腎長期存活率仍不理想，排斥反應(yīng)仍是影響移植腎長期存活的重要因素。巨噬細(xì)胞在T細(xì)胞介導(dǎo)的排斥反應(yīng)（T cell-mediated rejection，TCMR）、抗體介導(dǎo)的排斥反應(yīng)（antibody-mediated rejection，AMR）及慢性排斥反應(yīng)中均有明顯的浸潤[13-14]，然而目前關(guān)于巨噬細(xì)胞M1亞型在排斥反應(yīng)中的作用研究仍然不充分。因此本研究通過生物信息學(xué)的方法，分析微陣列數(shù)據(jù)集，鑒定了排斥反應(yīng)中和巨噬細(xì)胞M1亞型相關(guān)的基因，并基于6個M1-DEG構(gòu)建了預(yù)測移植物存活的風(fēng)險模型。

巨噬細(xì)胞具有高度可塑性，其在排斥反應(yīng)中發(fā)揮了多種作用，包括抗原加工和提呈、細(xì)胞毒性、吞噬作用及免疫反應(yīng)的調(diào)控等[7]。M1亞型作為促炎細(xì)胞的譜系，一方面可以分泌趨化因子以及細(xì)胞因子促進(jìn)炎癥[15]，另一方面可以促進(jìn)一氧化氮的產(chǎn)生引起組織損傷，最終影響移植腎的功能[16]。本研究通過WGCNA以及差異分析，在GSE36059數(shù)據(jù)集初步篩選出14個M1-DEG。隨后，在GSE21374的訓(xùn)練集中使用LASSO-Cox回歸篩選出TLR8、FCGR1B、BCL2A1、CTSS、GBP2及CARD16構(gòu)建風(fēng)險模型。TLR8是Toll樣受體家族中的重要成員，與病原體的識別以及先天免疫密切相關(guān)，可在腎移植術(shù)后誘導(dǎo)無菌性炎癥反應(yīng)并引起T細(xì)胞的募集，最終引起排斥反應(yīng)[17]。一項包括36例腎移植術(shù)后急性排斥反應(yīng)受者的研究，活組織檢查結(jié)果提示TLR8信使核糖核酸（messenger RNA，mRNA）表達(dá)上調(diào)[18]，與本研究發(fā)現(xiàn)一致。TLR8可能通過活化M1亞型以促進(jìn)其釋放炎癥因子、產(chǎn)生一氧化氮，從而損傷移植物[19]。

FCGR1B是Fcγ受體Ⅰ家族的成員，其通過結(jié)合免疫球蛋白γ的Fc區(qū)發(fā)揮作用。多中心的前瞻性研究發(fā)現(xiàn)腎移植受者血液中的FCGR1B的表達(dá)量與AMR的發(fā)生相關(guān)[20]，AMR受者外周血中FCGR1B的mRNA表達(dá)顯著上調(diào)。既往研究發(fā)現(xiàn)小鼠的巨噬細(xì)胞和樹突狀細(xì)胞高表達(dá)FCGR1[21]，因此FCGR1一方面可以與免疫球蛋白結(jié)合在體液免疫中發(fā)揮作用，另一方面可以通過介導(dǎo)非T細(xì)胞依賴的免疫復(fù)合物的攝取，聯(lián)系適應(yīng)性免疫和固有免疫。BCL2A1和細(xì)胞凋亡的進(jìn)程密切相關(guān)，在使用HLA預(yù)制抗體處理腎微血管內(nèi)皮細(xì)胞的研究中，加入IFN-γ和TNF-α組中BCL2A1的蛋白水平表達(dá)升高[22]。此外，在一項應(yīng)用支氣管肺泡灌洗液研究細(xì)胞因子風(fēng)暴的單細(xì)胞分析中，BCL2A1和巨噬細(xì)胞的炎癥反應(yīng)以及趨化作用相關(guān)[23]，進(jìn)一步表明該基因可能通過影響巨噬細(xì)胞的募集過程而影響移植物的存活。CTSS編碼的蛋白是肽酶C1家族的成員，其參與了抗原蛋白的降解呈遞給MHCⅡ類分子。一項小鼠心臟移植研究中使用了CTSS抑制劑，結(jié)果發(fā)現(xiàn)實驗組小鼠移植物存活時間更長，且沒有產(chǎn)生供體特異性抗體（donor specific antibody，DSA）[24]。同樣，在一項異基因小鼠脾臟移植研究中，移植脾臟組織CTSS表達(dá)量顯著上升，并且應(yīng)用CTSS抑制劑減少了T細(xì)胞的活化[25]。除了通過抗原提呈活化CD4+T細(xì)胞，CTSS也可以調(diào)控固有免疫誘導(dǎo)巨噬細(xì)胞釋放白細(xì)胞介素（interleukin，IL）-1β[26]。GBP2與IFN-γ信號傳導(dǎo)途徑有關(guān)，既往研究表明移植腎組織中該基因表達(dá)量可以用于診斷腎移植術(shù)后排斥反應(yīng)[27]。在一項包括46例肝移植受者的研究中，發(fā)生急性細(xì)胞性排斥反應(yīng)（acute cellular rejection，ACR）的受者外周血白細(xì)胞中GBP2 mRNA表達(dá)顯著上調(diào)，且多因素分析也表明GBP2高表達(dá)是ACR的獨(dú)立危險因素[28]。CARD16編碼的蛋白與半胱氨酸天冬氨酸蛋白酶（cysteinyl aspartate specific proteinase，Caspase）的激活以及隨后的IL-1β釋放相關(guān)。既往研究認(rèn)為CARD16是Caspase的抑制劑，抑制了后續(xù)的炎癥發(fā)展[29]。然而也有研究報道寡聚化的CARD16是Caspase的激活劑[30]。因此，CARD16可能通過影響巨噬細(xì)胞分泌IL-1β參與到移植物的損傷過程中?；?個M1-DEG的風(fēng)險模型對于預(yù)測移植腎的存活具有良好的表現(xiàn)，生存分析表明高、低風(fēng)險組的受者移植腎的存活存在明顯差異。因此，這些生物標(biāo)志物有助于在早期對腎移植受者進(jìn)行風(fēng)險分層，確定移植腎丟失高風(fēng)險受者，從而給予搶先干預(yù)，改善遠(yuǎn)期預(yù)后。

在高風(fēng)險組中，靜息及活化的CD4+記憶性T細(xì)胞、γδT細(xì)胞及巨噬細(xì)胞M1亞型浸潤增多。巨噬細(xì)胞通過將抗原提呈給CD4+T細(xì)胞使其活化，活化的CD4+T細(xì)胞通過產(chǎn)生IFN-γ和IL-2驅(qū)動TCMR，并產(chǎn)生IL-4、IL-5和IL-13參與AMR[31]。γδT細(xì)胞是固有免疫的重要組成部分，該細(xì)胞和移植后的排斥反應(yīng)及巨細(xì)胞病毒（cytomegalovirus，CMV）感染相關(guān)[32]。既往研究發(fā)現(xiàn)在腎移植術(shù)后排斥反應(yīng)期間移植腎組織中γδT細(xì)胞明顯浸潤，占總T細(xì)胞的10%以上[33]。在腎移植術(shù)后CMV感染后，循環(huán)中產(chǎn)生的CMV特異性γδT細(xì)胞有助于病毒的清除[34]。這些γδT細(xì)胞高表達(dá)CD16，因此其可能通過抗體依賴性細(xì)胞毒性參與由DSA促進(jìn)的同種異體移植損傷[35]。值得注意的是，在高風(fēng)險組的移植腎組織中M1亞型浸潤增多，而M0以及M2亞型減少，這揭示了巨噬細(xì)胞亞群平衡的失調(diào)是移植腎損傷的重要機(jī)制。GSEA分析進(jìn)一步揭示了在高風(fēng)險組中，免疫反應(yīng)、排斥反應(yīng)、細(xì)胞黏附分子及趨化因子信號通路明顯富集。最后，使用數(shù)據(jù)庫預(yù)測了調(diào)控預(yù)后基因的miRNA，這將為后續(xù)進(jìn)一步分析預(yù)后基因如何影響M1亞型的極化并造成移植物損傷提供依據(jù)。

綜上所述，本研究鑒定了排斥反應(yīng)中巨噬細(xì)胞M1亞型相關(guān)的基因，并構(gòu)建了風(fēng)險模型，對于預(yù)測移植腎存活具有良好的表現(xiàn)，為早期鑒定高風(fēng)險受者提供依據(jù)，但仍存在一些局限性。微陣列數(shù)據(jù)集無法闡明巨噬細(xì)胞極化過程中更加詳細(xì)的時間以及空間特征[6]，將巨噬細(xì)胞簡單地劃分為M1和M2兩個亞型無法滿足生物體內(nèi)不同組織以及病理狀態(tài)[36]。供者以及受者來源的巨噬細(xì)胞可能在排斥反應(yīng)中發(fā)揮了不同的作用，本研究無法區(qū)分巨噬細(xì)胞的來源，且需要后續(xù)的動物實驗進(jìn)一步驗證分子機(jī)制。