基于生物信息学和机器学习探究溃疡性结肠炎能量代谢关键基因及其潜在机制

李微, 李春梦, 刘青松, 张怡, 刘娅欣, 廖瑶玎

(成都中医药大学附属医院, 成都 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中的具体机制尚不明确,需要进一步的研究。

目前机器学习已经广泛应用于各种疾病的诊断、治疗以及预后,机器学习算法结合生物信息学分析为探寻疾病机制提供了新的思路和方法。现从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属于酰基辅酶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的发病机制和临床治疗提供新的思路和方向,并为基因功能验证、药物筛选和机制探索等基础实验指明了方向。