基于生物信息學(xué)方法篩選系統(tǒng)性紅斑狼瘡疾病進展中與先天免疫細胞相關(guān)的關(guān)鍵基因

左旭，李津

（天津醫(yī)科大學(xué)基礎(chǔ)醫(yī)學(xué)院生物信息學(xué)系，天津 300070）

系統(tǒng)性紅斑狼瘡（systemic lupus erythematosus，SLE）是一種最為常見的自身免疫性疾病。成年女性的發(fā)病率約為男性的8 倍，其發(fā)病累及全身多種器官，對患者的身心健康造成了巨大的危害[1]。SLE 的發(fā)病和進展是遺傳和環(huán)境因素共同作用的結(jié)果，并且受到雌激素水平的影響。有研究報道SLE 患者的一級親屬患病風(fēng)險比健康人群高20 倍[2]；環(huán)境因素，如紫外線和陽光照射會導(dǎo)致細胞凋亡增加，進而釋放富含自身抗原的細胞外囊泡，從而對SLE 的發(fā)病有一定的影響[3]；性別和激素影響也是SLE 的危險因素之一，已有研究發(fā)現(xiàn)，雌激素會通過刺激B 細胞增殖和隨后的自身抗體產(chǎn)生來促進體液免疫反應(yīng)。

目前，對SLE 發(fā)病機制的主要研究集中于致病性自身抗體及免疫復(fù)合物、T 細胞和自然殺傷（NK）細胞功能失調(diào)等方面[4]。總的來說，SLE 的發(fā)展主要是在各種病因的影響下導(dǎo)致致病性自身抗體的產(chǎn)生；另一方面部分SLE 患者由于CD8+T 細胞和NK細胞的功能失調(diào)，不能抑制CD4+T 細胞的持續(xù)作用，導(dǎo)致B 細胞持續(xù)活化產(chǎn)生自身抗體，繼而與細胞組織或循環(huán)中的抗原結(jié)合，形成免疫復(fù)合物通過沉積在各器官的組織中或直接引起組織和器官的損傷。先天性和適應(yīng)性免疫系統(tǒng)在SLE 的發(fā)病機制中都發(fā)揮作用[5－6]。異常的先天免疫反應(yīng)在SLE 的發(fā)病機制中起重要作用，其通過釋放炎性細胞因子導(dǎo)致組織損傷以及自身反應(yīng)性T 和B 細胞的異常激活，進而導(dǎo)致致病性自身抗體的產(chǎn)生和終末器官傷害。而SLE 患者中先天免疫細胞在狼瘡發(fā)病進展中的作用和機制還缺乏深入研究[7－8]。探索SLE 發(fā)生、發(fā)展的分子特征和機制，鑒定出與SLE 相關(guān)的遺傳學(xué)特征，為SLE 的有效預(yù)防、診斷和治療提供新的策略顯得非常重要。

微陣列數(shù)據(jù)（microarray data）以及RNA 高通量測序（RNA－sequencing，RNA－seq）數(shù)據(jù)已廣泛用于在基因組水平上探索和鑒定用于疾病早期診斷和預(yù)后的生物標志物[9]。有遺傳學(xué)研究表明，先天免疫過程在SLE 發(fā)病機制中也發(fā)揮重要作用[10]。例如，中性粒細胞溶解因子1（NCF1）和NCF2 基因中的錯義單核苷酸多態(tài)性（SNPs）與SLE 疾病相關(guān)[11]；單核細胞上FcγRI/CD64 的表達增加與SLE 的持續(xù)炎癥和腎炎相關(guān)[12]。但對于SLE 的先天免疫過程中的關(guān)鍵基因了解還不夠充分。因此在本研究中，基于微陣列數(shù)據(jù)以及RNA－Seq 數(shù)據(jù)，在轉(zhuǎn)錄水平分析出先天免疫細胞的異常與SLE 的發(fā)病高度相關(guān)，并通過加權(quán)共表達網(wǎng)絡(luò)分析（Weighted Gene Co－ex pression Network Analysis，WGCNA）構(gòu)建了共表達網(wǎng)絡(luò)，結(jié)合免疫細胞組成和豐度分析，篩選出與SLE先天免疫相關(guān)的樞紐基因，為SLE 的早期診斷及治療提供一定的幫助。

1 材料與方法

1.1 數(shù)據(jù)來源及預(yù)處理 從高通量基因表達數(shù)據(jù)庫（Gene Expression Omnibus data base，GEO）中共收集來自GSE50772、GSE81622 和GSE99967 的89例SLE 患者血液樣本和62 名健康對照血液樣本的基因表達微陣列數(shù)據(jù)，作為本研究的發(fā)現(xiàn)隊列。分別使用R（版本3．5．2）中Affy 包的rma（Robust Multiarray Average）行歸一化，并且對于有多個探針或探針組的基因取中位數(shù)的值。然后，采用KNearest Neighbors 的方法來處理3 個數(shù)據(jù)集中的缺失值。之后，首先使用R 軟件包的inSilicoMerging 來合并數(shù)據(jù)集，然后，使用sva 包的combat 功能去除批次效應(yīng)。選取GPL169791 Illumina HiSeq 2500（Homo sapiens）和GPL18573 Illumina NextSeq 500（Homo sapiens）平臺的GSE122459 數(shù)據(jù)集作為驗證隊列，其包括20 例SLE 患者血液樣本和6 名健康對照血液樣本。

1.2 免疫細胞組成和豐度分析 采用CIBERSORT和xCELL 方法對這89 例SLE 患者數(shù)據(jù)進行分析，CIBERSORT 采用反卷積算法通過具有22 種免疫細胞亞型的參考集（LM22）估計免疫細胞組成和豐度，xCELL 采用ssGSEA 算法通過計算樣本在每種細胞類型標記的富集分數(shù)并將其轉(zhuǎn)換為細胞類型分數(shù)，最后進行校正得到各細胞類型的xCELL 評分[13]。之后，篩選出先天免疫細胞并采用秩和檢驗來比較SLE 組和對照組之間的免疫細胞組成差異，使用R 中的ggplot2 包來進行可視化。

1.3 加權(quán)基因共表達網(wǎng)絡(luò)的構(gòu)建 通過利用89 例SLE 樣本的基因表達譜，剔除方差最小的前70%的基因，得到3 689 個基因，然后，使用R 中的WGCNA 包來構(gòu)建共表達網(wǎng)絡(luò)。首先，基于配對基因之間的Pearson 相關(guān)值，將單個轉(zhuǎn)錄本的表達水平轉(zhuǎn)換為相似矩陣。接下來，將相似度矩陣轉(zhuǎn)換為鄰接矩陣，計算公式為amn=|cmn|β（cmn=配對基因之間的Pearson 相關(guān)性；amn=配對基因之間的鄰接）。參數(shù)β可以提高基因間的強相關(guān)性，降低基因間的弱相關(guān)性。選擇β=9 的冪時，鄰接矩陣轉(zhuǎn)換為拓撲重疊矩陣（Topological overlap matrix，TOM）。用基于TOM不同性測量的平均連桿分層聚類對顯示相似表達譜的基因進行分類，這些基因模塊由簇樹的分支和不同顏色表示[14]。計算模塊特征基因與先天免疫細胞組成和豐度水平之間的相關(guān)性，通過Pearson 檢驗確定模塊的顯著性。

1.4 富集分析 使用R 軟件的clusterprofile 包進行GO（Gene Ontology）和KEGG（Kyoto Encyclopedia of Genes and Genomes，京都基因與基因組百科全書）途徑富集分析，以檢測與關(guān)鍵模塊基因相關(guān)的生物功能和潛在信號通路途徑。GO 分析確定生物過程（biological process，BP）、細胞組分（Cellular component，CC）、分子功能（molecular function，MF），P＜0．05 為篩選標準。KEGG 分析確定潛在信號通路途徑，F(xiàn)DR＜0．05 作為統(tǒng)計顯著性標準。

1.5 差異表達基因分析及關(guān)鍵基因的篩選 健康樣本和SLE 樣本間的差異表達基因（differentially expressed gene，DEG）是通過R 軟件包limma 來識別的，以|log2 FoldChange|≥1 和FDR＜0．05 作為有統(tǒng)計學(xué)意義的差異標準進行篩選。計算與基因的表達相關(guān)性以獲得基因顯著性（gene significance，GS），同時計算模塊特征向量與基因的表達相關(guān)性以獲得模塊成員關(guān)系（module membership，MM），選取了MM＞0．8，且GS＜0．6 的基因作為樞紐基因。此外，在關(guān)鍵模塊和DEG 中的基因之間重疊了基因以進一步確定關(guān)鍵基因。

1.6 基因集富集分析（gene set enrichment analysis，GSEA） 根據(jù)關(guān)鍵基因的表達水平的中位數(shù)將樣本分成高表達組和低表達組，以MSigDB 數(shù)據(jù)庫的c2．cp．kegg．v7．4．symbols．gmt 為參考基因集，用以評估相關(guān)途徑和分子機制，基于基因表達譜分組，設(shè)定最小基因集為5，最大基因集為5 000，經(jīng)1 000次重抽樣，以|NES|＞1，P＜0．05，F(xiàn)DR＜0．25 為明顯富集基因集的標準。

1.7 受試者工作特征（ROC）分析 利用R 軟件的pROC 包進行SLE 診斷中具有高靈敏度和特異性的關(guān)鍵基因的鑒定。并在GSE122459 數(shù)據(jù)集中計算了關(guān)鍵基因的表達與疾病中性粒細胞比例的斯皮爾曼相關(guān)系數(shù)。

2 結(jié)果

2.1 去批次效應(yīng)以及免疫細胞組成和豐度分析 將數(shù)據(jù)集進行整合，并進行去除批次效應(yīng)（圖1），然后通過CIBERSORT 和xCELL 進行免疫細胞組成和豐度分析，將先天免疫細胞的結(jié)果進行可視化。CIBERSORT 的分析發(fā)現(xiàn)SLE 組與對照組相比，靜息的NK 細胞和靜息的肥大細胞的比例減少，而單核細胞、巨噬細胞M0、活化的樹突狀細胞和中性粒細胞的比例增加（圖2A）；xCELL 的分析顯示，SLE組中同樣是NK 細胞和肥大細胞的比例較對照組少，而單核細胞、巨噬細胞、樹突狀細胞和中性粒細胞的比例更多（圖2B）。

圖1 去批次效應(yīng)UMAP 圖Fig 1 UMAP diagram before and after removing batch effect

圖2 免疫細胞組成和豐度分析Fig 2 Immune cell composition and abundance analysis

2．2 WGCNA 構(gòu)建共表達網(wǎng)絡(luò) 選取擬合指數(shù)R2為0．86 時的軟閾值β=9（圖3A），基于選擇的軟閾值構(gòu)建共表達網(wǎng)絡(luò)，共識別到13 個模塊，并通過層次聚類樹來展示所構(gòu)建的13 個模塊（圖3B）。

圖3 軟閾值的確定及共表達網(wǎng)絡(luò)構(gòu)建Fig 3 Determination of soft thresholds and co-expression network construction

2.3 免疫細胞相關(guān)模塊的確定及其KEGG 通路分析 將免疫細胞組成和豐度分析中均具有統(tǒng)計學(xué)意義的結(jié)果作為性狀與模塊內(nèi)的基因相關(guān)聯(lián)，結(jié)果顯示棕褐色模塊與靜息的NK 細胞（cor=0．80，P＜0．05）和綠色模塊與中性粒細胞（cor=0．85，P＜0．05）具有顯著相關(guān)性（圖4），故將綠色模塊和棕褐色模塊作為關(guān)鍵模塊進行后續(xù)分析。

圖4 模塊基因與性狀相關(guān)性圖Fig 4 Correlation between module genes and traits

2.4 富集分析 對棕褐色模塊的基因和綠色模塊的基因進行KEGG 富集分析，綠色模塊基因富集在白細胞介素（IL）－17 信號通路、核因子（NF）－κ B 信號轉(zhuǎn)導(dǎo)途徑、腫瘤壞死因子（TNF）信號轉(zhuǎn)導(dǎo)通路等（圖5A），棕褐色模塊基因富集在抗原處理和呈現(xiàn)、NK 細胞的細胞毒性等免疫系統(tǒng)相關(guān)的通路上（圖5C）。對棕褐色模塊和綠色模塊的基因進行GO 通路富集分析，綠色模塊基因主要參與免疫反應(yīng)、趨化作用和炎癥等生物途徑；在細胞組分（CC）中主要富集到分泌顆粒膜、分泌囊泡等；在分子功能（MF）中主要富集在細胞因子活性、信號受體結(jié)合等活動上（圖5B）。棕褐色模塊基因在生物過程（BP）中主要富集到免疫反應(yīng)、細胞殺傷等；CC 分析表明，這些基因產(chǎn)物的位置主要富集在細胞溶解顆粒；MF 富集顯示，它們的分子功能（MF）主要富集到MHC Ⅰ類受體活性等（圖5D）。

圖5 KEGG 和GO 富集分析Fig 5 KEGG and GO enrichment analysis

2.5 樞紐基因的篩選 計算綠色模塊和棕褐色模塊內(nèi)基因的GS 和MM 值并繪制散點圖，篩選了MM＞0．8、GS＞0．6 的基因作為樞紐基因，其中綠色模塊中有15 個樞紐基因（圖6A），棕褐色模塊中有8個樞紐基因（表1 和圖6）。

表1 綠色模塊和棕褐色模塊中的樞紐基因Tab 1 Hub genes in green module and tan module

圖6 關(guān)鍵模塊基因相關(guān)性散點圖Fig 6 Module eigengenes in the key module

2.6 差異表達基因分析及關(guān)鍵基因的篩選 對合并的數(shù)據(jù)集進行差異表達基因分析，得到上調(diào)差異基因45 個，下調(diào)差異基因5 個（圖7A）。將這些差異表達基因與篩選得到的樞紐基因整合，在綠色模塊中得到兩個關(guān)鍵基因，分別是CXCL1（C－X－C Motif Chemokine Ligand 1）和MME（Membrane Metalloendopeptidase）（圖7B）。

圖7 差異基因表達與關(guān)鍵基因篩選Fig 7 Differential gene expression and key gene screening

2.7 GSEA 分析 結(jié)果表明，CXCL1 富集到NOD樣受體信號通路（NES=1.97，P =0.002）、利什曼病的感染通路（NES=1.87，P=0.004）且均與其表達正相關(guān)，而硒胺酸的代謝通路（NES=-1.83，P=0.006）與其負相關(guān)（圖8A）。MME 富集到葉酸的代謝（NES=1.79，P=0.002）、Toll 樣受體信號通路（NES=1.82，P=0.002）和NOD 樣受體通路（NES=1.82，P=0.002）且均與MME 的表達呈正相關(guān)（圖8B）。

圖8 基因集富集分析Fig 8 Gene set enrichment analysis

2.8 GSE122459 數(shù)據(jù)集驗證 使用GSE122459 數(shù)據(jù)集對兩個關(guān)鍵基因CXCL1 與MME 進行驗證，在表達上，兩個基因在SLE 患者中表達量升高（圖9A），CXCL1 基因的AUC（Area under curve）值為0.86，MME 基因的AUC 值為0.81（圖9B）。對數(shù)據(jù)集進行免疫細胞比例分析后，計算其表達與中性粒細胞相關(guān)性，CXCL1 的表達與中性粒細胞的相關(guān)性為0.79（圖9C），MME 的表達與中性粒細胞的相關(guān)性為0.87（圖9D）。

圖9 CXCL1 和MME 在GSE122459 數(shù)據(jù)集中的驗證Fig 9 Validation of CXCL1 and MME in GSE122459 dataset

3 討論

SLE 具有高度復(fù)雜性和異質(zhì)性。患者免疫系統(tǒng)多組成部分出現(xiàn)明顯的功能障礙，并且呈現(xiàn)臨床異質(zhì)性，所以對SLE 的精準治療面臨著巨大的挑戰(zhàn)。關(guān)于SLE 發(fā)病過程中先天免疫系統(tǒng)的作用機制的報道較少，但最近有研究發(fā)現(xiàn)中性粒細胞鐵死亡對于人自身免疫性疾病有著重要作用，證實了先天免疫細胞異常會導(dǎo)致系統(tǒng)性自身免疫性疾病[15]。因此，開展對于SLE 發(fā)病過程中先天免疫的作用機制的研究具有重要意義。

本研究整合了GSE50772、GSE81622 和GSE99967數(shù)據(jù)集，進行了先天免疫細胞組成和豐度的全面評估，發(fā)現(xiàn)NK 細胞和肥大細胞在SLE 中的比例減少，而單核細胞、巨噬細胞、樹突狀細胞和中性粒細胞在SLE 中的比例較對照組增加，提示這些細胞類型在SLE 的發(fā)生和發(fā)展中發(fā)揮重要作用。之前的研究表明，NK 細胞的細胞毒性功能在SLE 患者中受損[16]；肥大細胞及其激活相關(guān)抗體參與了類風(fēng)濕性關(guān)節(jié)炎和多發(fā)性硬化癥等各種自身免疫性疾病的發(fā)生和發(fā)展[17]；低密度粒細胞（LDG，low-density granulocytes）即中性粒細胞的一種促炎性中性粒細胞亞群與狼瘡的特定臨床特征的存在和嚴重程度相關(guān)等[18]。

為了進一步鑒定與疾病相關(guān)的模板，對數(shù)據(jù)集進行了WGCNA 分析，共得到13 個基因共表達模塊，并鑒定了1 個與NK 細胞高度相關(guān)的模塊即棕褐色模塊和1 個與中性粒細胞高度相關(guān)的模塊即綠色模塊。與NK 細胞相關(guān)的模塊基因主要富集中在細胞溶解顆粒、細胞殺傷、抗原處理和呈現(xiàn)、NK細胞的細胞毒性等通路中；與中性粒細胞相關(guān)的模塊基因主要富集于細胞因子活性、IL-17 信號通路、NF-kappa B 信號傳導(dǎo)途徑等。有研究證實，對細胞因子的監(jiān)測能夠提高確定SLE 疾病活性的敏感性和特異性，如通過IL-18 能夠預(yù)測活動性腎臟SLE的風(fēng)險，而IL-6 和IL-8 能夠預(yù)測活動性非腎臟的風(fēng)險[19]。中性粒細胞在SLE 中發(fā)揮至關(guān)重要的作用。狼瘡相關(guān)免疫復(fù)合物在FcγRIIA 依賴但非TLR的響應(yīng)中激活人類中性粒細胞[20]，中性粒細胞外陷阱NET 中的mtDNA、抗mtDNA 抗體與SLE 中的PDC IFNα 發(fā)病機制之間相關(guān)聯(lián)[21]。

之后，結(jié)合DEG 分析選擇CXCL1 和MME 作為與中性粒細胞相關(guān)的關(guān)鍵基因。GSEA 分析發(fā)現(xiàn)，CXCL1 和MME 都有富集到NOD 樣受體信號通路。近期有研究報道NOD 樣受體信號通路及其通路相關(guān)基因在SLE 中的重要作用，NOD 樣受體含吡啶域蛋白3（NLRP3）炎癥體的失調(diào)在系統(tǒng)性紅斑狼瘡中起著重要作用，高水平的let-7f-5p 可以通過靶向NLRP3 來減輕SLE 炎癥[22]。此外，CXCL1 的表達還可能與利什曼病的感染通路正相關(guān)并與硒氨酸的代謝呈負相關(guān)，內(nèi)臟利什曼病（Visceral leishmaniasis，VL）與SLE 的癥狀有很強的相似性[23]；體內(nèi)的許多硒蛋白參與內(nèi)源性抗氧化防御系統(tǒng)，硒蛋白的水平低與各種疾病的發(fā)展密切相關(guān)[24]。MME 的表達還與葉酸的代謝和Toll 樣受體通路呈正相關(guān)。有研究表明，低水平的葉酸與高濃度的同型半胱氨酸相關(guān)，而高同型半胱氨酸血癥是SLE 患者亞臨床動脈粥樣硬化的獨立風(fēng)險因素之一[25-26]。網(wǎng)狀或凋亡中性粒細胞釋放的核酸通過病毒核酸特異性Toll 樣受體激活先天和適應(yīng)性免疫是狼瘡腎炎的發(fā)病機制之一[27]。

然后，基于GSE122459 數(shù)據(jù)集進行了ROC 分析以及基因的表達與中性粒細胞的相關(guān)性分析來評估基因與患者診斷之間的相關(guān)性。CXCL1（AUC=0.86）和MME（AUC=0.81）均在SLE 診斷中表現(xiàn)出較高的靈敏度和特異性，且其表達均與中性粒細胞有較高的相關(guān)性。CXCL1 屬于CXC 趨化因子家族的一員，可作為多種免疫細胞的趨化劑，是一種有效的中性粒細胞趨化劑和激活劑[28]，同時，CXCL1也能作為將活躍的白細胞招募到炎癥組織的重要原因之一。已有研究報道，CXCL1 表達水平的變化可能是SLE 活動的潛在標志[29]。MME 編碼膜金屬內(nèi)肽酶也稱中性內(nèi)肽酶，主要存在于腎臟、肺、大腦和肝臟細胞的血漿膜中。盡管目前沒有其在SLE 中的研究，但有研究表明在與SLE 相關(guān)的膜性腎病中，中性內(nèi)肽酶是重要的自身抗原之一，在調(diào)節(jié)炎癥方面至關(guān)重要[30]。因此，推測CXCL1 和MME 可能在SLE 的疾病進展中發(fā)揮重要作用。本研究存在一定的局限性，由于使用公共數(shù)據(jù)進行分析，對于SLE的臨床異質(zhì)性無法獲取全面的組織數(shù)據(jù)進行分析，從而選擇了外周血的數(shù)據(jù)進行分析。外周血是SLE免疫系統(tǒng)的主要途徑，外周血單核細胞（PBMC）是最先啟動對目標器官的自身免疫過程的免疫細胞，因此，PBMC 的基因表達特征能夠在一定程度上揭示目標器官中免疫細胞的分子特征。在本研究所揭示的關(guān)鍵基因還需要在更全面的組織樣本隊列中進一步驗證[31]。

綜上所述，本研究基于免疫細胞組成與WGCNA 識別了2 個與SLE 先天免疫相關(guān)的關(guān)鍵模塊，以及2 個與中性粒細胞相關(guān)的關(guān)鍵基因即CXCL1和MME，在SLE 發(fā)生、發(fā)展的先天免疫中可能發(fā)揮著重要作用，值得深入研究。