李微, 李春夢, 劉青松, 張怡, 劉婭欣, 廖瑤玎

(成都中醫(yī)藥大學(xué)附屬醫(yī)院, 成都 610075)

潰瘍性結(jié)腸炎(ulcerative colitis,UC)是一種累及結(jié)直腸黏膜及黏膜下層為主的慢性非特異性、非感染性、炎癥性腸道疾病[1]。其主要臨床表現(xiàn)為腹痛,腹瀉,黏液膿血便,遷延難愈,不僅影響患者的生活質(zhì)量,還增加了結(jié)腸癌等嚴(yán)重并發(fā)癥的風(fēng)險(xiǎn)。近年來,UC的發(fā)病率和患病率在全球范圍內(nèi)呈上升趨勢[2-3]。目前,UC是研究的熱點(diǎn),現(xiàn)有機(jī)制發(fā)現(xiàn)其跟遺傳、微生物、免疫力、環(huán)境等因素有關(guān),但其發(fā)病機(jī)制仍不清楚[4]。因此,探尋UC新的發(fā)病機(jī)制對于指導(dǎo)診斷和探尋新的治療途徑具有重要的意義。

能量代謝是機(jī)體物質(zhì)代謝過程中伴隨的能量釋放、轉(zhuǎn)移和利用過程,是維持人體機(jī)能活動的重要環(huán)節(jié)。腸上皮細(xì)胞是極性細(xì)胞,維持其極性依賴高能量,三磷酸腺苷(adenosine-triphosphate triphosphate,ATP)匱乏時(shí)會導(dǎo)致細(xì)胞骨架改變,細(xì)胞間緊密連接減少,腸黏膜通透性增加,從而導(dǎo)致UC的發(fā)生[5-6]。Roediger等用“能量缺乏”一詞來描述UC患者結(jié)腸細(xì)胞的代謝表現(xiàn)和病變黏膜的組織學(xué)變化,他們發(fā)現(xiàn)UC患者丁酸代謝明顯減少,尤其是結(jié)腸和直腸,因此提出丁酸代謝減少,腸黏膜屏障功能減退,從而引起炎癥的發(fā)生,導(dǎo)致UC的發(fā)生[7]。目前,有關(guān)潰瘍性結(jié)腸炎與能量代謝相關(guān)研究已有較大進(jìn)展,Sünderhauf等[8]將氧化磷酸化酶系調(diào)節(jié)的P32結(jié)腸低表達(dá)與UC常見的線粒體功能障礙、杯狀細(xì)胞分化缺陷和黏液屏障形成受損聯(lián)系起來,認(rèn)為黏膜中P32/gC1qR/HABP1的缺失導(dǎo)致潰瘍性結(jié)腸炎能量不足,并影響杯狀細(xì)胞的分化。研究表明,IL-6和TNF-α等促炎細(xì)胞因子在UC患者的能量代謝中發(fā)揮重要作用,治療后靜息能量消耗(measured resting energy expenditure,MREE)也會減少[9],通過減少清除損傷的線粒體和過多的自由基,以對抗UC中的氧化應(yīng)激損傷,保護(hù)上皮細(xì)胞,對UC有較好的治療作用[10]。但是,能量代謝在UC中的具體機(jī)制尚不明確,需要進(jìn)一步的研究。

目前機(jī)器學(xué)習(xí)已經(jīng)廣泛應(yīng)用于各種疾病的診斷、治療以及預(yù)后,機(jī)器學(xué)習(xí)算法結(jié)合生物信息學(xué)分析為探尋疾病機(jī)制提供了新的思路和方法?，F(xiàn)從GEO數(shù)據(jù)庫中下載UC數(shù)據(jù)集,篩選出其中與能量代謝相關(guān)的基因(differentially expressed energy-related genes,DEERGs),并且對這些基因進(jìn)行富集分析,再利用LASSO和SVM-RFE鑒定DEERGs中的關(guān)鍵基因,并進(jìn)行生物信息學(xué)分析,為闡述UC的機(jī)制提供一個(gè)新的方向。

1 材料和方法

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

從GEO數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/)中以“ulcerative colitis”為關(guān)鍵詞檢索UC相關(guān)數(shù)據(jù)集,下載了GSE87466和GSE75214兩個(gè)數(shù)據(jù)集。GSE87466包括87個(gè)UC樣本和21個(gè)正常組織樣本,GSE75214包括97個(gè)UC樣本和11個(gè)正常組織樣本,作為基因表達(dá)驗(yàn)證數(shù)據(jù)集。此外,還使用了從人類代謝相關(guān)途徑數(shù)據(jù)庫中下載(https://reactome.org/)的594個(gè)與能量代謝相關(guān)基因,此數(shù)據(jù)庫是代謝通路專用數(shù)據(jù)庫[11]。

1.2 能量代謝相關(guān)基因篩選

從GSE87466數(shù)據(jù)集中提取能量代謝基因的表達(dá)量,通過limma包對能量代謝基因進(jìn)行差異分析,閾值設(shè)置為|log2FC|>0.585,設(shè)置P<0.05篩選得到DEERGs,其中FC表示倍性變化(fold change)。同時(shí),使用R軟件包pheatmap和ggplot2軟件包繪制DEERGs的熱圖和火山圖。

1.3 基因本體和京都基因與基因組百科全書富集分析

為了揭示DEERGs的潛在生物功能,使用R軟件包Cluster Profiler進(jìn)行基因本體(gene ontology,GO)和京都基因與基因組百科全書(Kyoto encyclopedia of genes and genomes,KEGG)富集分析,設(shè)置P<0.05時(shí)具有統(tǒng)計(jì)學(xué)意義,即富集較顯著,并且將GO和KEGG分析的前10個(gè)結(jié)果進(jìn)行可視化繪制出GO和KEGG富集分析圖。

1.4 機(jī)器學(xué)習(xí)篩選關(guān)鍵基因

通過R軟件包glmnet、e107進(jìn)行最小絕對收縮和選擇算子(the least absolute shrinkage and selection operator,LASSO)算法和支持向量機(jī)-遞歸特征消除(support vector machine-recursive feature elimination,SVM-RFE)分析鑒定關(guān)鍵基因。用VennDiagram將二者獲得的基因取交集,最后所得到的交集即UC能量代謝的關(guān)鍵基因。對關(guān)鍵基因繪制接受者操作特征(receiver operating characteristics,ROC)曲線并計(jì)算曲線下面積(area under the curves,AUC),以檢測關(guān)鍵基因的診斷效能。

1.5 單基因GSEA和GSVA

為鑒定每個(gè)關(guān)鍵基因的潛在功能,使用R軟件包BiocManager、Org.Hs.eg.db從The Molecular Signatures Database(MSigDB,https://www.gsea-msigdb.org/gsea/msigdb)下載選定的參考基因組。用R軟件包c(diǎn)lusterProfiler對關(guān)鍵基因進(jìn)行GSEA-GO途徑分析,使用R軟件包GSVA對每個(gè)關(guān)鍵基因進(jìn)行GSVA-KEGG分析,對基因GSVA評分分?jǐn)?shù)進(jìn)行t檢驗(yàn)分析,并用R軟件包enrichplot、ggpubr繪制相應(yīng)的GSEA和GSVA直方圖。

1.6 免疫浸潤分析

使用CIBESORT算法獲得GSE87466數(shù)據(jù)集中22種免疫浸潤細(xì)胞的比例,并使用R軟件包tidyverse將5個(gè)關(guān)鍵基因?qū)?2種免疫滲透細(xì)胞的影響進(jìn)行相關(guān)性分析進(jìn)行可視化。

1.7 靶向藥物預(yù)測

從藥物-基因相互作用數(shù)據(jù)庫(DGIdb,https://dgidb.org)中搜索這5個(gè)關(guān)鍵基因的靶向藥物,并下載結(jié)果,將其導(dǎo)入Cytospace3.9.1軟件繪制出靶向藥物網(wǎng)絡(luò)圖譜。

1.8 構(gòu)建ceRNA網(wǎng)絡(luò)

使用R軟件包spongeScan、TargetScan和miRDB來預(yù)測5個(gè)關(guān)鍵基因結(jié)合的miRNAs,然后用R軟件包spongeScan找到與上述miRNAs結(jié)合的lncRNA存在的靶向關(guān)系,再導(dǎo)入Cytospace3.9.1軟件以構(gòu)建網(wǎng)絡(luò)圖。

1.9 關(guān)鍵基因表達(dá)驗(yàn)證

將關(guān)鍵基因在數(shù)據(jù)集GSE75214中進(jìn)行表達(dá)水平驗(yàn)證,并用R軟件包ggpubr對基因表達(dá)結(jié)果繪制差異箱線圖。

2 結(jié)果

2.1 差異分析結(jié)果

下載了594個(gè)能量代謝基因,差異分析后共篩選出32個(gè)DEERGs,包括9個(gè)上調(diào)基因和23個(gè)下調(diào)基因,使用火山圖(圖1)和熱圖(圖2)進(jìn)行可視化,并對這些DEERGs進(jìn)行了相關(guān)性分析(圖3)。

紅色代表上調(diào)基因;綠色代表下調(diào)基因圖1 DEERGs表達(dá)火山圖Fig.1 Volcano map of DEERGs expression

2.2 GO和KEGG富集分析

為了研究這些DEERGs的生物學(xué)功能,對其進(jìn)行GO和KEGG富集分析。GO富集分析(圖4)表明,DEERGs主要涉及外源性代謝過程、細(xì)胞對外源性刺激的反應(yīng)、胞內(nèi)葡萄糖醛酸化、尿酸代謝過程和脂肪酸代謝過程等通路有較為富集的結(jié)果。對于KEGG的富集分析結(jié)果(圖5)顯示,DEERGs與藥物代謝-細(xì)胞色素P450、視黃醇代謝、糖酵解/葡萄糖新生,抗壞血酸和醛代謝等代謝途徑富集結(jié)果較顯著。

圖4 GO富集分析Fig.4 GO enrichment analysis

圖5 KEGG富集分析Fig.5 KEGG enrichment analysis

2.3 關(guān)鍵基因的篩選

使用LASSO和SVM-RFE兩種機(jī)器學(xué)習(xí)鑒定了UC能量代謝相關(guān)關(guān)鍵基因,使用LASSO算法篩選出8個(gè)關(guān)鍵基因(圖6、圖7),然后使用SVM-RFE算法鑒定了8個(gè)關(guān)鍵基因(圖8、圖9)。最后,通過Venn圖譜取交集(圖10),篩選出5個(gè)(SLC16A1、ACSF2、CBR3、NR1H4、CHST11)作為UC能量代謝關(guān)鍵基因,構(gòu)建ROC曲線(圖11)并計(jì)算曲線下面積(AUC),發(fā)現(xiàn)這5個(gè)關(guān)鍵基因的AUC均大于0.9。

圖6 LASSO篩選繪制垂直線的最佳參數(shù)Fig.6 LASSO Filter and draw the best parameters for vertical lines

圖7 能量代謝相關(guān)基因的LASSO系數(shù)譜Fig.7 LASSO coefficient profiles of energy metabolism related genes

圖8 SVM交叉驗(yàn)證準(zhǔn)確性圖形Fig.8 SVM cross-validation accuracy graph

圖9 SVM交叉驗(yàn)證誤差圖形Fig.9 SVM cross-validation error graph

圖10 維恩圖Fig.10 Venn diagram

圖11 5個(gè)關(guān)鍵基因的ROC曲線Fig.11 ROC curves of the 5 key genes

2.4 基因富集分析

對這5個(gè)關(guān)鍵基因進(jìn)行了GSEA-GO功能富集分析(圖12～圖16),ACSF2富集調(diào)節(jié)細(xì)胞因子產(chǎn)生、炎癥反應(yīng)的正調(diào)節(jié)。CBR3富集功能為適應(yīng)性免疫反應(yīng)、粒細(xì)胞趨化性、淋巴細(xì)胞介導(dǎo)的免疫,CHST11和SLC16A1與B細(xì)胞受體信號通路功能較富集,NR1H4的各免疫途徑較為富集。

圖12 ACSF2的單基因功能富集分析Fig.12 Single gene functional enrichment analysis of ACSF2

圖13 CBR3的單基因功能富集分析Fig.13 Single gene functional enrichment analysis of CBR3

圖14 CHST11的單基因功能富集分析Fig.14 Single gene functional enrichment analysis of CHST11

圖15 NR1H4的單基因功能富集分析Fig.15 Single gene functional enrichment analysis of NR1H4

圖16 SLC16A1的單基因功能富集分析Fig.16 Single gene functional enrichment analysis of SLC16A1

2.5 基因變異富集分析(GSVA)

分析這5個(gè)關(guān)鍵基因高表達(dá)和低表達(dá)水平時(shí)的GSVA-KEGG(圖17～圖21)。ACSF2低表達(dá)時(shí),丁酸代謝、萜類骨架生物合成、丙酮酸代謝較富集;NR1H4低表達(dá)時(shí)外源性細(xì)胞色素P450代謝、視黃醇代謝、淀粉和蔗糖代謝較富集;SLC16A1低表達(dá)時(shí)候富集的信號通路是外源性細(xì)胞色素p450代謝、視黃醇代謝和淀粉和蔗糖代謝;CBR3高表達(dá)時(shí),纈氨酸-亮氨酸和異亮氨酸降解、丙酸代謝途徑較富集;CHST11高表達(dá)時(shí)丙酸代謝、丁酸代謝等途徑較富集。

圖17 ACSF2的單基因變異富集分析Fig.17 Single gene variation enrichment analysis of ACSF2

圖18 CBR3的單基因變異富集分析Fig.18 Single gene variation enrichment analysis of CBR3

圖19 CHST11的單基因變異富集分析Fig.19 Single gene variation enrichment analysis of CHST11

圖20 NR1H4的單基因變異富集分析Fig.20 Single gene variation enrichment analysis of NR1H4

圖21 SLC16A1的單基因變異富集分析Fig.21 Single gene variation enrichment analysis of SLC16A

2.6 免疫浸潤分析

免疫浸潤分析結(jié)果(圖22、圖23)表明UC與巨噬細(xì)胞M0、巨噬細(xì)胞M1和中性粒細(xì)胞呈正相關(guān),與巨噬細(xì)胞M2、肥大細(xì)胞負(fù)相關(guān)。將關(guān)鍵基因與免疫細(xì)胞進(jìn)行相關(guān)性分析發(fā)現(xiàn),巨噬細(xì)胞M0的表達(dá)與SLC16A和NR1H4呈負(fù)相關(guān),與CHST11和CBR3呈正相關(guān),中性粒細(xì)胞與CHST111呈正相關(guān),與CBR3和ACSF2呈負(fù)相關(guān),調(diào)節(jié)性T細(xì)胞的表達(dá)與CHST11呈負(fù)相關(guān),與ACSF2呈正相關(guān)。

圖22 UC中免疫細(xì)胞相關(guān)性Fig.22 Correlation of immune infiltrating cells in UC

2.7 靶向藥物預(yù)測

DGIdb數(shù)據(jù)庫未收錄ACSF2、SLC16A1、CHST11的藥物信息,共找到5種CBR3相關(guān)靶點(diǎn)藥物和30種NR1H4相關(guān)靶點(diǎn)藥物,將其導(dǎo)入Cytoscape3.9.1軟件中繪制圖譜后可視化結(jié)果(圖24)。

圖24 靶向藥物圖譜Fig.24 Targeted drug profiles

2.8 ceRNA網(wǎng)絡(luò)結(jié)果

構(gòu)建的ceRNA網(wǎng)絡(luò)結(jié)果見圖6D,該網(wǎng)絡(luò)含有265個(gè)節(jié)點(diǎn)(4個(gè)特征基因,111個(gè)miRNA和150個(gè)lncRNA)和297個(gè)邊,miR-182-5p和hsa-miR-590-3p、hsa-miR-149-3p對基因的調(diào)控作用較大(圖25)。

圖25 ceRNA網(wǎng)絡(luò)Fig.25 ceRNA network graph

2.9 數(shù)據(jù)集表達(dá)驗(yàn)證

最后,在數(shù)據(jù)GSE75214中驗(yàn)證了關(guān)鍵基因的表達(dá),結(jié)果表明,上述5個(gè)基因與數(shù)據(jù)集GSE87466中的表達(dá)一致。其中,ACSF2、NR11H4、SLC16A1在UC組織中的表達(dá)顯著降低。CBR3和CHST11的表達(dá)顯著高于正常組織(圖26)。P<0.01,表示差異具有統(tǒng)計(jì)學(xué)意義。

圖26 基因表達(dá)驗(yàn)證Fig.26 Gene expression validation

3 討論與結(jié)論

共篩選出了5個(gè)與UC中能量代謝相關(guān)的關(guān)鍵基因(SLC16A1、ACSF2、NR1H4、CHST11和CBR3),對其進(jìn)行GSVA和GSEA分析,結(jié)果表明5個(gè)關(guān)鍵基因主要與糖酵解/葡萄糖新生、丁酸代謝、丙酮酸代謝、淀粉和蔗糖代謝等代謝途徑以及B細(xì)胞受體信號通路、適應(yīng)性免疫反應(yīng)、細(xì)胞因子活性等炎癥調(diào)節(jié)途徑有著密切的關(guān)系。

SLC16A1,也稱為單羧酸鹽轉(zhuǎn)運(yùn)蛋白1(monocarboxylic acid transporter,MCT1),是一種人類單羧酸鹽轉(zhuǎn)運(yùn)體,通過人類結(jié)腸細(xì)胞的頂膜轉(zhuǎn)運(yùn)丁酸和乳酸[12]。當(dāng)SLC16A1的表達(dá)被RNA干擾抑制時(shí),丁酸誘導(dǎo)的細(xì)胞周期阻滯和分化也會減少[13]。人類結(jié)腸細(xì)胞的能量主要由葡萄糖代謝和氧化以及3種短鏈脂肪酸的乙酸鹽、丙酸鹽和丁酸鹽代謝提供,SLC16A1的下調(diào)與炎癥性結(jié)腸細(xì)胞對丁酸的利用受損有關(guān)[14]。

ACSF2屬于?；o酶A合成酶家族,它通過與輔酶A形成硫酯來催化脂肪酸代謝中的初始反應(yīng)[15]。Chen等[16]發(fā)現(xiàn)ACSF2可以作為腸上皮干細(xì)胞的脂肪酸氧化基因。同樣,Zhao等[17]發(fā)現(xiàn)ACSF2可以作為肝細(xì)胞癌的治療靶點(diǎn)和免疫相關(guān)生物標(biāo)志物。

CBR3是一種細(xì)胞質(zhì)和單體還原型輔酶Ⅱ 依賴性羰基還原酶,屬于短鏈脫氫酶超家族。Malatkova等[18]發(fā)現(xiàn),在炎癥環(huán)境的刺激下,CBR3的表達(dá)會被激活,是炎癥刺激的新靶點(diǎn)基因。

CHST11(碳酸氫鹽硫轉(zhuǎn)移酶11),也稱為軟骨素-4-硫轉(zhuǎn)移酶-1(C4ST-1),是參與軟骨素硫酸化的酶之一,通過催化硫酸鹽供體的硫酸基轉(zhuǎn)移到N-乙酰半乳糖的C-4位[19]。

NR1H4可以編碼一種激活法尼醇X受體(farnesoid X receptor,FXR)的轉(zhuǎn)錄因子,目前被認(rèn)為是一種膽汁酸受體和一種合成生物受體,參與整個(gè)膽汁酸代謝途徑,在維持腸道完整性方面發(fā)揮重要作用,與炎癥性腸病關(guān)系密切[20]。因此,能量代謝的關(guān)鍵基因可能通過脂肪酸代謝、葡萄糖代謝、脂質(zhì)代謝等各種代謝方式參與UC的發(fā)生和發(fā)展。

免疫浸潤分析顯示,這些基因與巨噬細(xì)胞、中性粒細(xì)胞和CD4+T細(xì)胞、Tregs細(xì)胞等免疫浸潤細(xì)胞密切相關(guān)。免疫細(xì)胞浸潤是UC生物學(xué)過程中不可或缺的因素;在UC患者的血液和結(jié)腸組織中發(fā)現(xiàn)大量活化的中性粒細(xì)胞,UC患者的巨噬細(xì)胞轉(zhuǎn)化為M1型,釋放IL-6、IL-12、TNF-α等促炎細(xì)胞因子以及活性氧,抑制細(xì)胞增殖并損傷周圍組織,而巨噬細(xì)胞M2通過釋放IL-10和TNF-β促進(jìn)細(xì)胞增殖和組織修復(fù)[21-23],此外,T細(xì)胞同樣也與UC的發(fā)生密切相關(guān),當(dāng)Treg/Th17平衡失調(diào)而導(dǎo)致免疫穩(wěn)態(tài)失衡時(shí)會引起UC發(fā)生[24]。因此,能量代謝關(guān)鍵基因還可能通過調(diào)控免疫浸潤的方式UC中發(fā)揮作用。

目前競爭內(nèi)源性RNA在UC的發(fā)病機(jī)制中得到重視,環(huán)狀RNA可以通過MiR-182-5p促進(jìn)腸細(xì)胞自噬來減輕UC患者的腸黏膜損傷[25]。miR-375-3p可被腸上皮細(xì)胞激活,并通過靶向激活UC中的STAT3信號傳導(dǎo)加劇炎癥反應(yīng)[26]。研究表明miR-133α參與人類結(jié)腸上皮細(xì)胞和實(shí)驗(yàn)性結(jié)腸炎的促炎信號傳導(dǎo)[27], miR-23a-3p和miR-27a-3p下調(diào)可以觸發(fā)結(jié)腸組織中線粒體代謝通路,破壞腸道黏膜屏障引起UC[28],因此,挖掘UC相關(guān)能量代謝的環(huán)狀RNA可作為診斷治療靶點(diǎn)新的切入點(diǎn)。

對能量代謝關(guān)鍵基因的潛在藥物進(jìn)行預(yù)測,發(fā)現(xiàn)可作用于NR1H4的藥物中,托非索作為一種FXR受體激動劑,已被證明可用于治療原發(fā)性膽管炎患者[29],而香膠甾酮可通過干擾NF-κB通路和炎癥因子來調(diào)節(jié)腸道炎癥,它為UC的治療提供了新的思路[30]。FXR激動劑奧貝膽酸可以通過抑制細(xì)菌移位來恢復(fù)腸道屏障并抑制腸道炎癥[31]。5種CBR3相關(guān)藥物用于治療各種腫瘤。其他還未被研究證實(shí)的藥物可能具有治療UC的潛力,值得進(jìn)一步的臨床驗(yàn)證。

近年來,機(jī)器學(xué)習(xí)算法結(jié)合生物信息學(xué)分析來研究復(fù)雜疾病的診斷預(yù)測和預(yù)后評估成為研究的熱點(diǎn)[32-33]。通過機(jī)器算法篩選出了5個(gè)與UC相關(guān)的能量代謝關(guān)鍵基因,這些基因可以通過脂肪酸代謝、葡萄糖代謝等代謝途徑和腸道免疫微環(huán)境的調(diào)節(jié)來影響UC的發(fā)生發(fā)展,為UC的發(fā)病機(jī)制和臨床治療提供新的思路和方向,并為基因功能驗(yàn)證、藥物篩選和機(jī)制探索等基礎(chǔ)實(shí)驗(yàn)指明了方向。