整合生物信息學(xué)分析痤瘡的關(guān)鍵基因與免疫浸潤水平

樊文龍，王紅心，陳東宇，楊曉雨，黃巧，胡敏，潘素躍，王樸，何玉清,2

1.廣東醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系，廣東醫(yī)科大學(xué)醫(yī)學(xué)系統(tǒng)生物學(xué)研究所，廣東 東莞 523808；2.東莞市寮步醫(yī)院，廣東 東莞 523400

痤瘡作為一種生活中常見的慢性炎癥性皮膚病，在不同人群中均有發(fā)生，好發(fā)于青春期[1]，在面部、胸背部較為常見，表現(xiàn)為粉刺、丘疹、膿皰、結(jié)節(jié)、囊腫甚至瘢痕等不同程度的損害。2019年，全球痤瘡發(fā)病人數(shù)約為1.17億例，1990—2019年期間，全球痤瘡的整體患病率和疾病負(fù)擔(dān)均逐年升高[2]。痤瘡可能會(huì)影響患者生活質(zhì)量[3-5]。痤瘡的發(fā)生受多種因素影響，主要包括遺傳、雄激素刺激皮脂腺異常角化、微生物影響，炎癥的病理免疫反應(yīng)等[6]。雄激素/雄激素樣受體、Toll樣受體、表皮生長因子受體、過氧化物酶體增殖物激活受體、成纖維細(xì)胞生長因子受體和FOXO1轉(zhuǎn)錄因子等均參與尋常痤瘡的發(fā)病過程[7]，但其發(fā)病機(jī)制仍不完全明晰。近年來，利用生物信息學(xué)的方法進(jìn)行數(shù)據(jù)挖掘已然成為一項(xiàng)探究發(fā)病機(jī)制的研究手段。因此，本研究通過生物信息學(xué)方法對(duì)從GEO數(shù)據(jù)庫獲取的數(shù)據(jù)進(jìn)行分析，探究與痤瘡發(fā)病相關(guān)的關(guān)鍵基因，為后續(xù)相關(guān)研究提供理論參考。

1 材料與方法

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

從GEO數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/)中下載三個(gè)數(shù)據(jù)集：GSE108110(GPL570平臺(tái))、GSE53795(GPL570平臺(tái))和GSE6475(GPL571平臺(tái))；其中GSE108110更新于2020年3月；GSE53795更新于2019年3月；GSE6475更新于2018年12月。這三個(gè)數(shù)據(jù)集的樣本均為痤瘡患者皮損組織和正常組織，其中GSE108110皮損組織和正常組織各18例，GSE53795皮損組織和正常組織各12例，GSE6475皮損組織和正常組織各6例，共納入72例樣品組織。

1.2 差異表達(dá)基因(differentially expressed genes, DEGs)的提取

利用Rstudio對(duì)下載好的數(shù)據(jù)進(jìn)行處理，并安裝“impute” R包對(duì)缺失值采用最近鄰居法(K-nearest neighbor, KNN)進(jìn)行填充，安裝 “l(fā)imma” R包對(duì)處理好的基因進(jìn)行分析，在校正后P<0.05，|logFC|>1的標(biāo)準(zhǔn)下，篩選出痤瘡的DEGs，繪制韋恩圖篩選出三個(gè)數(shù)據(jù)集的共同DEGs。使用在線分析軟件平臺(tái)SangerBox(http://vip.sangerbox.com/)，去除三個(gè)數(shù)據(jù)集的批次效應(yīng)并合并數(shù)據(jù)，用于后續(xù)分析，該平臺(tái)整集100多種常用的分析方法，涵蓋統(tǒng)計(jì)、分析、可視化等工具，提供了大量廣泛使用的生物信息學(xué)分析工具[8]。

1.3 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis, WGCNA )[9]

根據(jù)標(biāo)準(zhǔn)化后的基因表達(dá)矩陣，篩選出皮損組和對(duì)照組有差異的基因(P<0.05)，利用SangerBox構(gòu)建WGCNA，計(jì)算軟閾值，進(jìn)行基因聚類和模塊劃分(最小模塊基因數(shù)取30，敏感性取3，模塊合并閾值取0.25)，評(píng)價(jià)基因模塊和痤瘡表型的關(guān)聯(lián)性，篩選出正相關(guān)和負(fù)相關(guān)最高的模塊。

1.4 基因功能和通路富集分析

對(duì)兩基因模塊中的基因和共同DEGs取交集，使用SangerBox對(duì)交集基因進(jìn)行基因本體論(gene ontology, GO)功能注釋，同時(shí)對(duì)交集基因進(jìn)行京都基因與基因組百科全書(kyoto encyclopedia of genes and genomes, KEGG)通路富集分析，在P<0.05的條件下，篩選顯著的功能和通路。

1.5 蛋白質(zhì)互作網(wǎng)絡(luò)(protein-protein interaction network, PPI network)的構(gòu)建與分析

使用STRING在線工具(https://string-db.org/)分析蛋白質(zhì)-蛋白質(zhì)相互作用，并使用Cytoscape(version:3.8.2)進(jìn)行可視化，構(gòu)建PPI網(wǎng)絡(luò)[10-11]；使用CytoHubba插件中的4種常用算法(MCC、Degree、Betweenness和Closeness)分別篩選出PPI網(wǎng)絡(luò)中前20的關(guān)鍵基因，并取它們的交集。

1.6 免疫細(xì)胞浸潤分析

基于SangerBox在線分析軟件，通過CIBERSORT算法進(jìn)一步研究痤瘡樣本中22種免疫細(xì)胞浸潤的差異[12]，評(píng)估22種免疫細(xì)胞浸潤的比例，分析關(guān)鍵基因與免疫細(xì)胞浸潤的相關(guān)性。

2 結(jié)果

2.1 差異表達(dá)基因提取

數(shù)據(jù)處理后，在校正后P<0.05，|logFC|>1的條件下，GSE108110、GSE53795和GSE6475數(shù)據(jù)集中分別篩出547、1 106和284個(gè)DEGs，其中GSE108110上調(diào)基因472個(gè)，下調(diào)基因75個(gè)；GSE53795上調(diào)基因981個(gè)，下調(diào)基因125個(gè)；GSE6475上調(diào)基因261個(gè)，下調(diào)基因23個(gè)。篩選三個(gè)數(shù)據(jù)集的共同DEGs，共篩選出154個(gè)DEGs，其中上調(diào)基因152個(gè)，下調(diào)基因2個(gè)(表1，圖1)。

2.2 WGCNA

WGCNA分析結(jié)果顯示，軟閾值β為24，signed R2為0.85，平均連接度為3.41，網(wǎng)絡(luò)符合無標(biāo)度網(wǎng)絡(luò)的分布。在最小模塊基因數(shù)為30、敏感性為3、模塊合并閾值為0.25的標(biāo)準(zhǔn)下，基因集被分為了4個(gè)共表達(dá)模塊，其中l(wèi)ightcyan模塊(4 419個(gè)基因)與樣本表型顯著正相關(guān)(r=0.92,P<0.001)，green模塊(228個(gè)基因)與樣本表型顯著負(fù)相關(guān)(r=-0.65,P<0.001，圖2)。

2.3 基因功能和通路富集

對(duì)WGCNA中正、負(fù)相關(guān)性最高的兩模塊進(jìn)行功能富集分析，GO富集結(jié)果顯示，在生物過程方面，主要集中在免疫系統(tǒng)過程、免疫反應(yīng)、細(xì)胞激活等；在細(xì)胞成分方面，基因主要集中于質(zhì)膜部分、囊泡、細(xì)胞內(nèi)囊泡等；在分子功能部分，基因主要與信號(hào)受體結(jié)合、含有蛋白質(zhì)的復(fù)合體結(jié)合、催化活性等有關(guān)(圖3A)；KEGG結(jié)果顯示，基因主要集中在細(xì)胞因子-細(xì)胞因子受體的相互作用、趨化因子信號(hào)傳導(dǎo)、NF-κB信號(hào)傳導(dǎo)等途徑(圖3B)。對(duì)這154個(gè)DEGs進(jìn)行功能富集分析，GO富集結(jié)果顯示，在生物過程方面，主要集中在免疫系統(tǒng)過程、免疫反應(yīng)、細(xì)胞激活等；在細(xì)胞成分方面，基因主要集中于分泌顆粒、分泌囊泡、細(xì)胞質(zhì)囊泡等；在分子功能部分，基因主要與信號(hào)受體結(jié)合、細(xì)胞因子活性、趨化因子受體結(jié)合等有關(guān)(圖3C)；KEGG結(jié)果顯示，基因主要集中在病毒蛋白與細(xì)胞因子和細(xì)胞因子受體的相互作用、細(xì)胞因子-細(xì)胞因子受體的相互作用、趨化因子信號(hào)傳導(dǎo)途徑、NF-κB信號(hào)傳導(dǎo)途徑、IL-17信號(hào)傳導(dǎo)等途徑(圖3D)。將WGCNA兩模塊基因與DEGs取交集，結(jié)果顯示，交集基因?yàn)?54個(gè)(圖3E)。

表1 共同差異表達(dá)基因Table 1 Differentially expressed genes

圖1 差異表達(dá)基因 1A、1B、1C：三個(gè)數(shù)據(jù)集的火山圖，紅色代表上調(diào)基因，綠色代表下調(diào)基因，黑色代表沒有顯著性差異的基因；1D：三個(gè)數(shù)據(jù)集上調(diào)基因的韋恩圖；1E：三個(gè)數(shù)據(jù)集下調(diào)基因的韋恩圖

圖2 各模塊與痤瘡的相關(guān)性Figure 2 Correlation of each module with acne.

2.4 蛋白質(zhì)互作網(wǎng)絡(luò)

使用Cytoscape構(gòu)建PPI網(wǎng)絡(luò)，得到139個(gè)節(jié)點(diǎn)，1 597條邊，提取連接度排名前10的節(jié)點(diǎn)，分別是IL-6、CXCL8、IL-1B、PTPRC、CD86、MMP9、ITGB2、CCR5、CCL2、TLR1(圖4A)。調(diào)用CytoHubba插件，采用4種算法篩選關(guān)鍵基因，提取4種算法的前20關(guān)鍵基因取交集，最終交集到2個(gè)關(guān)鍵基因，分別為CXCL8、CCR5(圖4B)。

圖3 基因功能和通路富集 3A、3B:兩模塊的GO富集和KEGG富集;3C、3D:DEGs的GO富集和KEGG富集;3E:兩模塊與DEGs的交集

圖4 PPI網(wǎng)絡(luò)(4A)及四種算法關(guān)鍵基因的交集(4B)

2.5 免疫細(xì)胞浸潤分析

使用CIBERSORT算法分析皮損組織與正常組織之間的差異，圖5顯示了兩組之間22 種免疫細(xì)胞類型的一般分布，與正常組織相比，皮損組織有更多成熟B細(xì)胞、靜息的CD4記憶T細(xì)胞、活化的CD4記憶T細(xì)胞、單核細(xì)胞、巨噬細(xì)胞M0、活化的樹突狀細(xì)胞、激活的肥大細(xì)胞和中性粒細(xì)胞，而記憶B細(xì)胞、漿細(xì)胞、調(diào)節(jié)性T細(xì)胞、靜息的樹突狀細(xì)胞、靜息的肥大細(xì)胞相對(duì)更少。

圖5 22種免疫細(xì)胞在痤瘡皮損組織和正常組織中浸潤的分布

3 討論

本研究采用生物信息學(xué)方法，對(duì)數(shù)據(jù)集GSE108110、GSE53795和GSE6475進(jìn)行數(shù)據(jù)挖掘，以便進(jìn)一步了解痤瘡的發(fā)病機(jī)制。通過WGCNA篩選重要基因模塊后，與差異表達(dá)基因取交集，最終得到了154個(gè)與痤瘡表型相關(guān)的差異表達(dá)基因，其中上調(diào)基因數(shù)量遠(yuǎn)高于下調(diào)基因，可能原因是痤瘡作為一種慢性炎癥性皮膚病，痤瘡發(fā)生后炎癥相關(guān)基因表達(dá)，同時(shí)本研究篩選差異基因設(shè)置的標(biāo)準(zhǔn)較為嚴(yán)格，導(dǎo)致某些下調(diào)基因被過濾掉。

進(jìn)行蛋白質(zhì)互作分析，采用4種常用的篩選關(guān)鍵基因的方法取交集，最終篩選出2個(gè)hub基因，CXCL8和CCR5。CXCL8 是C-X-C基序趨化因子配體，即白細(xì)胞介素-8(IL-8)，由上皮細(xì)胞和巨噬細(xì)胞等分泌。研究表明，CXCL8是腫瘤發(fā)生的關(guān)鍵因素，在腫瘤微環(huán)境中參與信號(hào)轉(zhuǎn)導(dǎo)，促進(jìn)腫瘤的進(jìn)展[13]，同時(shí)，CXCL8也參與多種皮膚病的發(fā)生相關(guān)，動(dòng)物實(shí)驗(yàn)表明，特應(yīng)性皮炎、酒渣鼻、銀屑病的患者中，CXCL8均有顯著增加[14-16]。另外，既往研究表明[17-18]，IL-8也參與痤瘡的發(fā)生發(fā)展，P.acnes通過激活Toll樣受體2(TLR2)能夠分泌IL-8，在SZ95皮脂細(xì)胞中，采用痤瘡丙酸桿菌無細(xì)胞提取物也能夠激活人NF-κB和p38 MAPK通路，并通過TLR2依賴的信號(hào)途徑上調(diào)IL-8的分泌，導(dǎo)致炎癥的發(fā)生，這與本研究中CXCL8作為上調(diào)基因的結(jié)果一致。CCR5為C-C基序趨化因子受體，在細(xì)胞趨化中發(fā)揮著重要作用，CCR5與其配體結(jié)合后，能夠促進(jìn)炎癥反應(yīng)以及誘導(dǎo)免疫反應(yīng)中不同T細(xì)胞亞群的黏附和遷移[19]。痤瘡是一種雄激素依賴性、胰島素樣生長因子IGF-1驅(qū)動(dòng)的皮脂腺毛囊疾病[20]，而相關(guān)文獻(xiàn)報(bào)道，CCR5過度表達(dá)與高雄激素血癥和胰島素抵抗有關(guān)[21]，CCR5的缺乏能夠?qū)е轮窘M織巨噬細(xì)胞中M2顯性表型的改變，有助于減輕肥胖誘導(dǎo)的胰島素抵抗[22-23]。叉頭盒轉(zhuǎn)錄因子O1(FOXO1)作為脂肪生成的負(fù)調(diào)控因子，能夠減少痤瘡的發(fā)生[24]，而痤瘡患者體內(nèi)胰島素和胰島素樣生長因子-1的水平升高會(huì)增加Akt的活性，從而通過磷酸化降低FOXO1的核活性[25]，增加痤瘡的發(fā)生，本研究也顯示CCR5上調(diào)，因此，推斷CCR5可能通過炎癥反應(yīng)和胰島素/胰島素生長因子-1信號(hào)調(diào)節(jié)參與痤瘡的發(fā)生。另外，本文KEGG富集結(jié)果顯示，篩選的兩個(gè)關(guān)鍵基因均富集于病毒蛋白與細(xì)胞因子和細(xì)胞因子受體的相互作用、細(xì)胞因子-細(xì)胞因子受體的相互作用、趨化因子信號(hào)傳導(dǎo)途徑，這兩個(gè)關(guān)鍵基因可能通過這些通路參與痤瘡的發(fā)生，而本文尚未進(jìn)行實(shí)驗(yàn)對(duì)篩選出的基因進(jìn)行驗(yàn)證，有待進(jìn)一步實(shí)驗(yàn)證實(shí)。

既往文獻(xiàn)表明，IL-17和Th17細(xì)胞與痤瘡發(fā)病有關(guān)，IL-17可作為痤瘡發(fā)病機(jī)制的標(biāo)志物，而且可能是痤瘡嚴(yán)重程度和瘢痕形成的潛在預(yù)后預(yù)測因子，Th17通路被激活后，可在疾病中發(fā)揮著關(guān)鍵作用[26-28]。Th17細(xì)胞屬于CD4+T淋巴細(xì)胞亞群，CD4+T細(xì)胞在轉(zhuǎn)化生長因子TGF-β與IL-21的共同誘導(dǎo)下分化為Th17細(xì)胞，而Th17細(xì)胞能夠特異性分泌細(xì)胞因子IL-17，IL-17可刺激角質(zhì)形成細(xì)胞產(chǎn)生促炎細(xì)胞因子、趨化因子和基質(zhì)金屬蛋白酶，引起炎癥反應(yīng)和組織破壞[29]。本研究在KEGG分析中富集到IL-17信號(hào)傳導(dǎo)通路，在一定程度上表明IL-17在痤瘡發(fā)病機(jī)制中占據(jù)重要地位，也間接表明CD4+T細(xì)胞的重要性，正如相關(guān)文獻(xiàn)所報(bào)道，在痤瘡初期，患者皮損及血清的CD4+T細(xì)胞較正常人明顯增多，經(jīng)過治療后，CD4+T細(xì)胞減少，痤瘡炎癥病變明顯改善[30-31]。本研究利用基因芯片數(shù)據(jù)集對(duì)22種免疫細(xì)胞浸潤分析顯示，相對(duì)于正常對(duì)照組，皮損組中有更多的中性粒細(xì)胞、活化的肥大細(xì)胞、活化的樹突狀細(xì)胞、活化的CD4記憶T細(xì)胞等，而調(diào)節(jié)性T細(xì)胞、靜息的樹突狀細(xì)胞、靜息的肥大細(xì)胞等相對(duì)更少，免疫細(xì)胞浸潤在痤瘡發(fā)生發(fā)展中的重要性不言而喻，但這些細(xì)胞具體的機(jī)制還需要進(jìn)一步驗(yàn)證。

綜上所述，本研究通過WGCNA和“l(fā)imma” R包篩選出了154個(gè)差異表達(dá)基因，主要集中在病毒蛋白與細(xì)胞因子和細(xì)胞因子受體的相互作用、細(xì)胞因子-細(xì)胞因子受體的相互作用、趨化因子信號(hào)傳導(dǎo)、NF-κB信號(hào)傳導(dǎo)、IL-17信號(hào)傳導(dǎo)等通路。CXCL8、CCR5這2個(gè)hub基因可能成為痤瘡發(fā)病的關(guān)鍵基因，免疫細(xì)胞浸潤水平的變化可能是痤瘡發(fā)病的重要原因，在痤瘡的發(fā)病中發(fā)揮著重要作用。這些結(jié)果有助于進(jìn)一步闡述痤瘡的發(fā)病機(jī)制，為后續(xù)的相關(guān)研究提供一定的參考。