摘要
基于生物信息学探讨炎症性肠病（Inflammatory bowel disease，IBD）与肌少症潜在的生物联系，为IBD合并肌少症的诊断和治疗提供新思路。从GEO数据库提取肌少症和IBD的数据集，筛选出共表达基因，并进行功能富集分析、蛋白质-蛋白质相互作用网络(PPI)、转录因子(TF)-基因相互作用网络和TF-miRNA协同调节网络的构建，最后通过DSigDB数据库预测药物分子。本研究共筛选到304个共表达基因，通过功能富集分析和PPI网络分析，得到3个关键基因（ICAM1、CXCL12、JUN），基于TF-基因相互作用网络和TF-miRNA相互作用网络分析，得出共21个miRNAs和17个TF共同参与肌少症和IBD的相互作用。在药物预测方面的分析结果表明，芹菜素、姜黄素等药物可能是治疗IBD合并肌少症的潜在药物分子。
Abstract
To explore the potential biological relationship between inflammatory bowel disease (IBD) and sarcopenia based on bioinformatics, and to provide a new idea for the diagnosis and treatment of IBD combined with sarcopenia. The data sets of sarcopenia and IBD were extracted from GEO database, co-expressed genes were screened, and functional enrichment analysis, protein-protein interaction network (PPI), transcription factor (TF) -gene interaction network and TF-miRNA co-regulatory network were constructed. Finally, drug molecules were predicted by DSigDB database. In this study, a total of 304 co-expressed genes were screened, and three key genes (ICAM1, CXCL12, JUN) were obtained through functional enrichment analysis and PPI network analysis. Based on TF-gene interaction network and TF-miRNA interaction network analysis, A total of 21 miRNAs and 17 TFS were involved in the interaction between sarcopenia and IBD. The results of drug prediction showed that apigenin, curcumin and other drugs may be potential drug molecules for the treatment of IBD combined with sarcopenia. By exploring the biological relationship between IBD and sarcopenia, we can provide new research ideas for the pathogenesis and further diagnosis and treatment of IBD with sarcopenia.
Keywords
炎症性肠病（Inflammatory bowel disease，IBD）是一种病因不明，免疫异常介导的慢性肠道炎症，常伴有复发性，包括溃疡性结肠炎和克罗恩病。IBD患者多伴有营养不良表现，主要决定因素是热量摄入减少、吸收不良、基础能量消耗增加、脂质和碳水化合物代谢发生显著改变[1]。然而，肌少症在临床上的发病率并不低，研究表明，与健康组对比，IBD患者骨骼肌明显减少[2-3], 患者的发生率在36.7%~65%之间[3]，严重影响患者的生活质量。肌少症是一种以骨骼肌质量、力量和功能进行性下降为特征的复杂综合征，该疾病现在已成为全球威胁。目前全球约有5 000万人患有肌肉减少症，预计到2050年，将高达5亿之多[4]。数据显示，65岁及以上人群中有1/3患有骨骼肌减少症(发病率为14%~33%)，80岁及以上老年人患病率高达50%~60%[4]。骨骼肌减少症原因复杂，不仅包括环境因素和遗传倾向，也包括环境和遗传的相互作用[5-6]。此外，有关遗传研究表明，肌肉质量和肌肉力量的遗传风险分别在32%~67%和46%~76%。而IBD的病因也是多因素的，包括大量的遗传背景[7-8]。然而，肌少症作为IBD中一种常见的并发症，其可能与IBD慢性炎症、营养缺陷和长期使用糖皮质激素在内的因素对肌肉和骨骼都有直接和间接的不良影响有关[9]，但确切发病机制尚无广泛研究。目前，治疗方案多为适当的营养和体育锻炼。因此，本研究将通过生物信息学技术，基于GEO数据库挖掘IBD与肌少症潜在的生物学途径和关系，同时预测IBD合并肌少症的潜在药物治疗靶点，为IBD合并肌少症的诊断和治疗提供新的方法和途径。
1 资料与方法
1.1 资料来源
数据库：GEO数据库（http://www.ncbi.nlm.nih.gov/geo），基因本体数据库（gene ontology, GO, http://geneontology.org/）, 京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG，https://www.genome.jp/kegg/），STRING数据库（http://www.string-db.org/），NetworkAnalyst在线平台（https://www.networkanalyst.ca/NetworkAnalyst/home.xhtml），JASPAR数据库（https://jaspar.genereg.net/），Enrichr平台（https://maayanlab.cloud/Enrichr/），DSigDB数据库（https://maayanlab.cloud/Enrichr/enrich#）。
1.2 方法
1.2.1 GEO数据库进行差异基因和GO/KEGG功能富集分析
在GEO数据库中，分别输入 “inflammatory bowel disease” 与 “sarcopenia” 进行检索，筛选出契合本研究的数据集。采用R语言（3.6.3版本），运用“limma”包中的normalizeBetweenArrays函数，对原始数据实施归一化处理，从而得到标准化的基因表达数据集。利用“ggplot2”包，针对GSE1428和GSE36807数据集进行分析，识别其中显著差异表达的基因。其筛选标准设定为|log2(FC)|>0.5且P<0.05，将满足此标准的基因确定为各数据集的差异表达基因，并通过Venn图展示。针对从Venn图中获取的共表达差异基因（DEGs），进一步开展GO和KEGG数据库分析，以挖掘其潜在生物学功能与相关通路。
1.2.2 基于STRING数据库构建DEGs蛋白-蛋白相互作用网络(PPI)和关键基因鉴定
在STRING 11.0数据库界面输入从Venn图中获得的DEGs，对DEGs构建差异表达基因的PPI网络，并通过Cytoscape软件将分析得到的PPI网络进行可视化分析。高度互连的关键基因由Cytoscape软件的插件cytoHubba鉴定。
1.2.3 NetworkAnalyst平台鉴定关键基因相关的转录因子和miRNA
在NetworkAnalyst平台界面输入关键基因，检索转录因子-基因和转录因子-miRNA相互作用网络，使用数据库默认分析方法，可视化共同调控网络。
1.2.4 基于DSigDB数据库预测关键基因的药物分子
首先在Enrichr平台界面输入关键基因，然后在Diseases/Drugs页面找到DSigDB数据库，检索与关键基因相关的预测药物分子。
1.3 统计学方法
采用数据库默认的统计学方法，使用R(3.6.3版本) 的“limma”包、“ggplot2”包、“clusterProfiler”包对数据进行差异表达和功能富集分析，P<0.05具有统计学意义。使用Cytoscape软件对DEGs进行可视化网络分析。
2 结果
2.1 DEGs的识别
GSE1428数据集包含10例年轻人和12例老年人的股外侧肌活检样本，共获得2 999个DEGs，其中有2 963个上调基因和36个下调基因符合筛选标准（图1(a)）。GSE36807数据集包含13个克罗恩病（Crohn’s disease，CD）样本、15个溃疡性结肠炎（Ulcerative colitis，UC）样本和7个健康对照样本，共获得符合筛选标准的1 530个DEGs，其中上调基因有739个，下调基因有791个（图1(b)）。使用R比较GES1428数据集和GSE36807数据集相互作用基因，共筛选到304个基因。使用韦恩图显示出肌少症和IBD的DEGs重叠信息（图1(c)）。
图
1
(a) GSE1428数据集火山图分析；(b) GSE36807数据集火山图分析；(c) GES1428数据集和GSE36807数据集Venn图分析
Fig.
1
(a) Volcano map analysis of GSE1428 data set; (b) Volcano map analysis of GSE36807 data set ; (c) Venn graph analysis of GES1428 dataset and GSE36807 dataset
2.2 GO/KEGG富集分析
基因功能富集分析使用R软件的clusterProfiler包进行分析。GO分析结果表明：常见的DEGs主要参与小分子代谢过程的调节、细胞外基质组织、细胞外结构组织、对氧气水平的反应和肽基酪氨酸磷酸化（图2(a)）；在分子功能组中，常见的DEGs主要参与蛋白酪氨酸激酶活性、跨膜受体蛋白激酶活性、受体配体活性、跨膜受体蛋白酪氨酸激酶活性和生长因子结合（图2(b)）；细胞成分研究表明，常见的DEGs主要参与含胶原蛋白的细胞外基质、细胞器外膜、外膜、线粒体外膜和富含ficolin-1颗粒（图2(c)）。基于KEGG的通路结果分析表明，常见的DEGs富集于阿米巴病、细胞因子-细胞因子受体相互作用、病毒蛋白与细胞因子和细胞因子受体的相互作用、Axon制导和蛋白酶体（图2(d)），各富集分析结果详见表1。
图
2
(a) GO富集分析生物过程结果分析图；(b) GO富集分析分子功能结果分析图；(c) GO富集分析细胞成分结果分析图；(d) KEGG富集分析结果图
Fig.
2
(a) Analysis of biological process results of GO enrichment analysis;(b) Analysis of molecular function results of GO enrichment analysis; (c) Analysis diagram of GO enrichment analysis of cell components; (d) Map of KEGG enrichment analysis results
表
1
GO/KEGG富集分析前五个信息列表
Table
1
List of top five information for GO/KEGG enrichment analysis
2.3 PPI网络分析和中枢基因识别
为了研究DEGs相对应的蛋白质之间的相互作用，使用STRING数据库（https://cn.string-db.org）构建了304个DEGs的PPI网络。使用Cytoscape软件可视化整个PPI网络（图3A）。PPI网络由241个节点和1 218条边组成。Cytoscape中的MCODE插件用于构建功能模块。结果表明，得分最高的模块（得分=7）由7个基因和21个边组成（图3B）。将Cytoscape中CytoHubba插件的10种拓扑方法用于筛选前10个中枢基因，其中MNC算法结果如图3C所示。在所有拓扑方法中发现3个共同基因，即ICAM1、CXCL12、JUN（表2）。
2.4 TF-gene相互作用网络
TF-gene相互作用网络是使用NetworkAnalyst网络工具生成的。对3个关键DEGs（ICAM1、CXCL12、JUN）确定TF-gene（图4）。该网络包含24个节点和25条边。ICAM1、CXCL12和JUN分别由10、9和6个TF-gene调控（转录因子靶点来源于JASPAR TF结合位点谱数据库）。
2.5 TF-miRNA协同调节网络
基于3个关键DEGs（ICAM1、CXCL12、JUN）的主要TF-miRNA协同调节网络通过NetworkAnalyst生成。创建的网络包括2 956个节点和2 993条边。通过筛选相互作用网络大于或等于2的TF和miRNA，进行可视化网络分析，该网络包含21个miRNAs和17个TF基因（图5）。
图
3
(a) DEGs的PPI网络图；(b) MCODE插件得分最高的模块分析图；(c) MNC算法分析结果图
Fig.
3
(a) PPI network of DEGs; (b) Module analysis diagram of MCODE plug-in with the highest score;(c) Figure of the results of MNC algorithm analysis
表
2
DEGs的10种拓扑算法分析结果
Table
2
Results of 10 topological algorithms analysis of DEGs
2.6 候选药物的鉴定
从Enrichr平台上的DSigDB数据库中检测药物分子，从候选药物中提取综合得分排在前15位的药物分子（表3）。这些药物可能是治疗肌少症和IBD的常见药物分子。结果表明，芹菜素 CTD 00007429和2,2'，4,4'，5,5'-六氯联苯 CTD 00000731是综合得分最高的两个药物分子。
图
4
关键DEGs（ICAM1、CXCL12、JUN）的TF-gene相互作用网络
Fig.
4
TF-gene interaction network of key DEGs (ICAM1, CXCL12, JUN)
图
5
关键DEGs（ICAM1、CXCL12、JUN）的TF-miRNA协同调节网络
Fig.
5
TF-miRNA co-regulatory network of key DEGs (ICAM1, CXCL12, JUN)
表
3
候选药物预测分子分析结果
Table
3
Results of predictive molecular analysis of candidate drugs
3 讨论
肌少症是IBD常见的并发症之一，在疾病活动、营养方面和糖皮质激素等方面，可通过多种机制影响IBD的肌肉骨骼的预后。既往的一项荟萃分析显示，约有42%的IBD患者合并肌少症，且肌少症在UC和CD患者中的发病率略有不同，其发病率分别约为37%和52%[10]。同样有研究表明与健康同龄人相比，IBD的肌肉质量和功能较低[11]。这些肌肉缺陷可能导致疾病的整体活动性降低，并与IBD治疗的反应性[12-13]呈负相关，因此，关注IBD合并肌少症的潜在机制，寻找其作用靶点，对于IBD疾病的有效治疗和提升患者的生活质量均有着重要的影响。本研究基于生物信息学分析技术，初步挖掘IBD合并肌少症的潜在相关性，旨在为其进一步研究发病机制及寻找诊治新靶点提供新的思路与方法。
在本研究共筛选到304个DEGs来分析肌少症和IBD之间的联系。基于GO功能富集结果显示，在生物功能过程中，部分DEGs显著参与小分子代谢过程的调节、细胞外基质组织、细胞外结构组织、对氧气水平的反应和肽基酪氨酸磷酸化，其中部分DEGs已在肌少症研究领域具有初步研究结果，如NF-kB-inducing激酶表达诱导E3泛素连接酶MuRF1和atrogin1 mRNA表达显著增加，进而会促进小鼠肌肉萎缩[14]。肌肉生长抑制素抑制细胞增殖基因的转录和促生长信号(如胰岛素和胰岛素样生长因子-1)的合成代谢通路，抑制蛋白质合成，从而强烈影响骨骼肌生长发育[15]。KEGG通路结果包括阿米巴病、细胞因子-细胞因子受体相互作用、病毒蛋白与细胞因子和细胞因子受体的相互作用、Axon制导和蛋白酶体。多种细胞因子-细胞因子-受体的相互作用在肌少症和IBD中得到证实[16-19]，对IBD患者骨骼肌形成具有负面影响，因此，细胞因子-细胞因子-受体的调控对IBD合并肌少症的发病机制有着潜在的重要作用。此外，丙酮酸脱氢酶激酶4的已知转录因子叉头盒蛋白O1通过转录调控MAFbx也可促进肌肉萎缩过程中的蛋白质降解，从而促进肌肉萎缩[20]。
PPI网络分析是本研究重要组成部分，有助于检测到与肌少症和IBD共同机制相关的潜在枢纽基因。本研究结果显示，3个关键基因分别是ICAM-1、CXCL12和JUN基因。有研究表明[21]ICAM-1修饰的间充质干细胞（MSCs）通过促进MSC归巢靶器官和免疫器官，显著减轻IBD小鼠的炎症损伤，是IBD治疗中新的潜在靶点，而其同时在肌少症中的表达，提示有可能通过营养干预改善肌少症而降低ICAM-1的表达[22]。CXCL12主要由骨髓、胸腺、淋巴结等组织中的基质细胞产生，参与炎症反应、肿瘤转移等病理过程中白细胞的迁移调控。G-蛋白偶联受体CXCR4与其天然配体CXCL12之间的相互作用与IBD患者所经历的炎症有关，阻断这种相互作用可能会减少IBD患者的炎症状态[23]。同时，有研究表明，CXCL12对骨代谢稳态调节产生影响，使得肌线粒体功能发生障碍，从而导致骨量的减少[24]。JUN是编码一种与病毒蛋白高度相似的蛋白，与特定的靶DNA序列直接相互作用，调节基因的表达。
本研究中基于TF-gene相互作用网络结果显示，在PPI分析中的关键基因（即ICAM1、CXCL12和JUN）分别有10,9和6个TF调控。更为重要的是，FOXC1同时调控ICAM1和JUN；TFAP2A同时调控JUN和CXCL12；NFKB1和RELA各调控了ICAM1和CXCL12两个基因。这就提示共同调控的基因有可能是IBD合并肌少症的关键分子，是疾病进展的关键环节。同时，TF-miRNA共调控网络分析通过可视化的TF-miRNA共调控网络显示了关键DEGs, TF和miRNA之间的关系，结果发现21个miRNAs和17个TF基因与ICAM1、CXCL12和JUN之间存在大于或等于2的相互调节网络。因此，这些关键基因（ICAM1、CXCL12、JUN）, TFs和miRNAs在IBD合并肌少症的发病机制存在着潜在相关性，为进一步研究其发病机制提供重要参考依据。
药物预测同样是本研究的重要部分。以DSigDB数据库为基础，检索关键DEGs（ICAM1、CXCL12、JUN）相关的药物化合物。在所有的预测药物中，我们的研究突出了前15名的重要药物。更为关键的是，关键DEGs（ICAM1、CXCL12、JUN）在前12种药物预测分子中均有提示。研究表明[25,26]芹菜素可以有效降低小鼠炎症细胞因子TNF-α、IL-1β、IL-6、MCP-1、CSF-1和COX-2的水平，可预防小鼠的慢性溃疡性结肠炎的发生。此外，芹菜素-7-O-葡萄糖苷可作为治疗肌少症的补充剂或补救途径[27]。此外，姜黄素类似物C66通过抑制DSS诱导的结肠组织中JNK与NF-κB信号通路的激活，可显著降低疾病活动指数和体重增加，可能成为治疗IBD的潜在候选药物[28]。同时，在小鼠实验中，姜黄素可维持老年再生比目鱼肌肌纤维成熟的成年水平和增加老年后肢肌肉中分离的MyoD阳性卫星细胞的百分比，影响卫星细胞。因此，姜黄素治疗通过改善卫星细胞来成功预防肌少症和肌少症的发展[29]。此外多项研究表明，姜黄素可通过抗细胞凋亡、抗炎和抗氧化介导的与肌肉降解和蛋白质合成相关的基因表达调节来预防肌肉疾病[30-32]。由此可见，这些预测药物分子在治疗肌少症和IBD中起着重要作用，需要进一步进行临床试验验证，并为今后的IBD合并肌少症的药物研发和治疗提供参考依据。
然而，尽管有这些初步的研究和理论基础，但具体的病理生理机制尚不清楚，需要更多的基础研究和临床研究来深入探讨IBD和肌少症之间的关系。特别是需要大规模的流行病学研究和病例对照研究来确认它们之间的关联性和可能的共同病因。
4 结论
本研究借助生物信息学综合分析手段，揭示了IBD与肌少症之间的潜在关系，得出以下结论：1）IBD与肌少症之间存在众多共同的关键基因和致病途径，为深入探索两者发病机制提供理论基础。2）研究基于关键基因预测出芹菜素、姜黄素等药物靶点，在缓解IBD炎症和抑制肌肉降解方面有潜在作用，为IBD合并肌少症的药物研发和治疗提供了新思路。
