基于生物信息學和機器學習探究潰瘍性結腸炎能量代謝關鍵基因及其潛在機制

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

(成都中醫藥大學附屬醫院, 成都 610075)

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

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

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

1 材料和方法

1.1 數據來源

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

1.2 能量代謝相關基因篩選

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

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

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

1.4 機器學習篩選關鍵基因

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

1.5 單基因GSEA和GSVA

為鑒定每個關鍵基因的潛在功能,使用R軟件包BiocManager、Org.Hs.eg.db從The Molecular Signatures Database(MSigDB,https://www.gsea-msigdb.org/gsea/msigdb)下載選定的參考基因組。用R軟件包clusterProfiler對關鍵基因進行GSEA-GO途徑分析,使用R軟件包GSVA對每個關鍵基因進行GSVA-KEGG分析,對基因GSVA評分分數進行t檢驗分析,并用R軟件包enrichplot、ggpubr繪制相應的GSEA和GSVA直方圖。

1.6 免疫浸潤分析

使用CIBESORT算法獲得GSE87466數據集中22種免疫浸潤細胞的比例,并使用R軟件包tidyverse將5個關鍵基因對22種免疫滲透細胞的影響進行相關性分析進行可視化。

1.7 靶向藥物預測

從藥物-基因相互作用數據庫(DGIdb,https://dgidb.org)中搜索這5個關鍵基因的靶向藥物,并下載結果,將其導入Cytospace3.9.1軟件繪制出靶向藥物網絡圖譜。

1.8 構建ceRNA網絡

使用R軟件包spongeScan、TargetScan和miRDB來預測5個關鍵基因結合的miRNAs,然后用R軟件包spongeScan找到與上述miRNAs結合的lncRNA存在的靶向關系,再導入Cytospace3.9.1軟件以構建網絡圖。

1.9 關鍵基因表達驗證

將關鍵基因在數據集GSE75214中進行表達水平驗證,并用R軟件包ggpubr對基因表達結果繪制差異箱線圖。

2 結果

2.1 差異分析結果

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

紅色代表上調基因;綠色代表下調基因圖1 DEERGs表達火山圖Fig.1 Volcano map of DEERGs expression

2.2 GO和KEGG富集分析

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

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

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

2.3 關鍵基因的篩選

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

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

圖7 能量代謝相關基因的LASSO系數譜Fig.7 LASSO coefficient profiles of energy metabolism related genes

圖8 SVM交叉驗證準確性圖形Fig.8 SVM cross-validation accuracy graph

圖9 SVM交叉驗證誤差圖形Fig.9 SVM cross-validation error graph

圖10 維恩圖Fig.10 Venn diagram

圖11 5個關鍵基因的ROC曲線Fig.11 ROC curves of the 5 key genes

2.4 基因富集分析

對這5個關鍵基因進行了GSEA-GO功能富集分析(圖12～圖16),ACSF2富集調節細胞因子產生、炎癥反應的正調節。CBR3富集功能為適應性免疫反應、粒細胞趨化性、淋巴細胞介導的免疫,CHST11和SLC16A1與B細胞受體信號通路功能較富集,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個關鍵基因高表達和低表達水平時的GSVA-KEGG(圖17～圖21)。ACSF2低表達時,丁酸代謝、萜類骨架生物合成、丙酮酸代謝較富集;NR1H4低表達時外源性細胞色素P450代謝、視黃醇代謝、淀粉和蔗糖代謝較富集;SLC16A1低表達時候富集的信號通路是外源性細胞色素p450代謝、視黃醇代謝和淀粉和蔗糖代謝;CBR3高表達時,纈氨酸-亮氨酸和異亮氨酸降解、丙酸代謝途徑較富集;CHST11高表達時丙酸代謝、丁酸代謝等途徑較富集。

圖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 免疫浸潤分析

免疫浸潤分析結果(圖22、圖23)表明UC與巨噬細胞M0、巨噬細胞M1和中性粒細胞呈正相關,與巨噬細胞M2、肥大細胞負相關。將關鍵基因與免疫細胞進行相關性分析發現,巨噬細胞M0的表達與SLC16A和NR1H4呈負相關,與CHST11和CBR3呈正相關,中性粒細胞與CHST111呈正相關,與CBR3和ACSF2呈負相關,調節性T細胞的表達與CHST11呈負相關,與ACSF2呈正相關。

圖22 UC中免疫細胞相關性Fig.22 Correlation of immune infiltrating cells in UC

2.7 靶向藥物預測

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

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

2.8 ceRNA網絡結果

構建的ceRNA網絡結果見圖6D,該網絡含有265個節點(4個特征基因,111個miRNA和150個lncRNA)和297個邊,miR-182-5p和hsa-miR-590-3p、hsa-miR-149-3p對基因的調控作用較大(圖25)。

圖25 ceRNA網絡Fig.25 ceRNA network graph

2.9 數據集表達驗證

最后,在數據GSE75214中驗證了關鍵基因的表達,結果表明,上述5個基因與數據集GSE87466中的表達一致。其中,ACSF2、NR11H4、SLC16A1在UC組織中的表達顯著降低。CBR3和CHST11的表達顯著高于正常組織(圖26)。P<0.01,表示差異具有統計學意義。

圖26 基因表達驗證Fig.26 Gene expression validation

3 討論與結論

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

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

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

CBR3是一種細胞質和單體還原型輔酶Ⅱ 依賴性羰基還原酶,屬于短鏈脫氫酶超家族。Malatkova等[18]發現,在炎癥環境的刺激下,CBR3的表達會被激活,是炎癥刺激的新靶點基因。

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

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

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

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

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

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