基于生物信息学分析溃疡性结肠炎的差异表达基因和相关微小RNA

吴胜男, 蒋环宇, 陈浩然, 王欣瑶, 吴佳辉, 王璐琦

(山西医科大学, 1. 管理学院, 2. 基础医学院, 山西 太原, 030000)

溃疡性结肠炎(UC)是一种慢性非特异性炎症,其特点是结直肠黏膜连续性、弥漫性发炎。UC作为炎症性肠病(IBD)的一个重要类型,主要表现包括腹泻、血便、腹痛和里急后重等症状,对患者的心理、经济状况和生活质量造成了严重的影响[1-3]。UC在工业化国家更为常见,特别是北美和西欧国家[4]。在IBD发病率迅速上升的亚太地区国家中,中国的发病率最高, UC是最常见的类型[5]。

目前UC的发病机制尚不明确,临床很难根治。UC的发病主要是由于肠道黏膜屏障功能的损伤、黏膜免疫失衡以及肠道菌群紊乱,同时参与发病的因素还有遗传易感性和环境因素[6]。免疫反应在UC的发生过程中扮演了重要角色,丧失免疫耐受会增加各种免疫炎症因子,刺激抗原特异性效应物的增殖,从而导致局部和全身炎症[7]。肠道损伤时,固有免疫细胞中性粒细胞能够迅速聚集到受损部位。中性粒细胞起到限制外源性微生物入侵的关键作用,并具备重要的免疫调节和杀菌功能,但过度的免疫反应却可能加剧UC早期肠道炎症和黏膜损伤[8]。相关研究[9]通过流式细胞术检测发现,与健康个体相比,UC患者的外周血T细胞表现出细胞激活标志物(如人类白细胞抗原-DR和β1整合素)的增加,而幼稚T细胞特征性抗原L选择素(CD62L)的表达下降,提示在UC患者中T细胞的活化程度增加。微小RNA(miRNA)可参与调节各种生理过程[10-11]。超过60%的人类蛋白质编码基因含有至少1个保守的miRNA结合位点[12]。miRNA与多种疾病(包括UC)具有显著相关性。

本研究对GEO数据库中的健康个体与UC患者的芯片原始数据进行加权基因共表达网络分析和差异分析,从而获取差异表达基因,并对其进行富集分析,根据关键的基因预测与DEG相关的潜在miRNA。

1 材料与方法

1.1 数据获取

通过NCBI的GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)以“ulcerative colitis” “ulcerative colitis clinic”为关键词进行检索,限制Entry type为“Series”、Top Organisms为“Homo sapiens”, 筛选之后获得GSE66407、GSE65114基因芯片数据作为研究对象,所在平台分别为GPL19833以及GPL16686。GSE66407中包含368个样本,通过筛选疾病以及是否具有表型数据,最终纳入的数据为对照组99个样本, UC组为156个样本; GSE65114中包含28个样本,其中对照组为12个, UC组为16个,全部纳入研究。

1.2 数据预处理及筛选差异表达基因

通过bioDBnet(https://biodbnet-abcc.ncifcrf.gov/)对GSE66407数据进行id转化,用于后续差异分析; 通过R软件中的dplyr包对GSE65114数据进行id转化,最终得到二者相对应的基因表达矩阵文件。对于GSE65114的基因表达矩阵文件,采用R软件中的limma包进行差异基因筛选。差异基因的筛选标准为|logFC|>1以及经过调整后P值<0.05。

1.3 基因集富集分析(GSEA)

本研究使用GSEA软件(version 4.2.3)对GSE65114数据进行分组。为评估相关的途径和分子机制,从MolecularSignatures Database(https://www.gsea-msigdb.org/gsea/msigdb)下载了c2.cp.kegg.v2023.1.Hs.symbols.gmt子集合。设置基因集的大小范围,最小基因集为15, 最大基因集为500, 并进行1 000次重复抽样。设定P值<0.05和FDR<0.25的阈值。通过基因表达谱数据和表型数据进行分析,得出了相关结果。

1.4 构建加权基因共表达网络(WGCNA)

当基因在生物学上显示高度相关时,使用WGCNA将其聚类到不同的模块中。此外, WGCNA还能够探索模块与外部性状之间的关联性,进而有助于识别候选的生物标志物和潜在的治疗靶点[13]。本研究的目的在于整合GSE66407的基因表达矩阵文件和临床表型文件,并利用R软件中提供的WGCNA包来构建WGCNA。本研究首先对数据进行聚类,以便检测异常值,并确定适当的软阈值来建立一个无尺度的网络。然后,采用层次聚类树和拓扑重叠矩阵的方法来鉴定和检测基因模块。为了确定各个模块与表型特征的相关性,本研究计算了皮尔逊相关系数,并找出与表型特征相关性最强的基因模块。

1.5 交集基因的功能富集分析

首先将差异基因与基因模块中的基因导入venny 2.1.0中进行交集分析,以获取交集基因。然后使用R软件中的ClusterProfiler包对这些交集基因进行功能富集分析,包括基因本体论(GO)和京都基因与基因组百科全书(KEGG)通路分析。GO功能分析包括生物学过程(BP)、细胞组分(CC)以及分子功能(MF)。设定P值<0.05作为筛选标准,用于确定显著富集的基因集条目。

1.6 蛋白互作网络的构建

为了揭示UC相关蛋白质的作用及其相互关系,将交集基因导入在线数据库STRING (https://cn.string-db.org/)中,以构建蛋白质相互作用(PPI)网络。在构建过程中,设定confidence>0.4作为筛选标准。随后,利用Cytoscape软件对构建好的PPI网络进行可视化操作。通过使用CytoNCA插件中Degree值的排序,选出排名前10位的基因作为枢纽基因(Hub基因)。

1.7 与Hub基因相关的miRNA

使用NetworkAnalyst 3.0(https://www.networkanalyst.ca/)在线平台[14], 将得到的10个Hub基因映射到各自的miRNA, 用来识别基因调控网络中miRNA-基因之间的相互作用。最后,通过Cytoscape软件对miRNA-基因网络图进行可视化。

2 结 果

2.1 差异基因的筛选结果

基于GSE65114的基因表达数据,将对照组和UC组进行比较。为了筛选差异基因并进行可视化,使用R软件中的ggplot2包。筛选标准为: |logFC|>1且P<0.05。最终,经过数据处理,得到了277个差异表达基因,其中有200个基因上调,77个基因下调。详细的基因热图和火山图见图1。在热图中,红色代表基因上调,绿色代表基因下调,而黑色则表示中间表达水平。基因表达越显著,颜色越深。

2.2 GSEA分析结果

通过对GSE65114数据中所有基因的表达矩阵进行GSEA分析,结果显示富集通路主要涉及神经活性配体受体相互作用通路、利什曼原虫感染、朊病毒病害以及心电图受体相互作用通路等,见图2。

2.3 基因共表达网络

本研究使用R软件中的WGCNA包对GSE66407的基因表达数据进行分析,首先利用goodSamplesGenes方法检测缺失值,并利用平均值计算基因的平均表达量,设置过滤标准为FPKM>0.5。通过WGCNA算法计算合适加权系数β(软阈值),得到β为11时构建的网络更符合无尺度特征(图3A)。接着通过层次聚类法,使用动态剪切树进行模块的识别和合并,最终获得16个共表达模块(图3B)。然后随机选取1 000个基因做TOM热图(图3C)发现了明显的颜色聚集,用来描述基因之间的拓扑重叠矩阵,确定每个模块独立存在。

2.4 模块与性状关联分析

本研究将样本临床特征中的年龄、炎症、组织纳入表型数据,将每个模块与样本的性状特征进行关联分析,发现brown模块中的基因与UC样本中发炎结肠的相关性最强(cor=0.69,P<0.01)(图4A)。

对统计棕色模块中基因的基因与性状关系(GS)和基因与模块关系(MM)进行分析,可以得出棕色模块内基因与显著性关系的散点图(图4C), 由散点图中可发现,基因的GS与MM在棕色模型中基本上都是呈现线性排布,存在着相关性与一致性,说明该模块与UC炎症的发生相关性最大,该模块中可能存在调控UC炎症反应的关键基因,因此之后将对棕色模块中的基因进一步分析。

2.5 交集基因的富集分析结果

利用GEO数据库中的GSE66407和GSE65114基因表达数据集进行分析,将获得的277个差异基因与基因模块中的2 027个基因导入venny 2.1.0中取交集,获得123个交集基因(图5A), 其中上调基因114个,下调基因9个。利用R软件中的ClusterProfiler包对所获得的交集基因进行GO功能和KEGG通路富集分析,设定P<0.05为富集条目的筛选标准。

GO富集结果(图5B)显示,与BP相关的条目主要富集在白细胞迁移、正向调节对外部刺激的反应、白细胞趋化性、髓系白细胞迁移、单核细胞迁移、中性粒细胞趋化性等;与CC相关的条目主要富集在胶原蛋白、细胞外基质分泌颗粒腔、细胞质囊泡腔、等离子体的囊泡腔外侧、细胞膜顶端部分等;与MF相关的条目主要富集在结构成分趋化因子活性、肝素结合、趋化因子受体结合、Toll样受体结合、细胞因子活性等。KEGG富集分析结果(图5C)显示,前12条通路主要富集在细胞因子受体相互作用、磷脂酰肌醇-3激酶/蛋白激酶B(PI3K-AKT)信号通路、趋化因子信号通路、核转录因子kappa B(NF-κB)信号通路、Toll样受体信号传导通路、白细胞介素-17(IL-17)信号通路等相关过程。

2.6 PPI网络及Hub基因的筛选结果

将交集基因导入String在线数据库,从中获取PPI网络(见图6A)。然后,利用Cytoscape软件对所获得的PPI网络进行可视化操作(见图6B), 使用CytoNCA插件计算节点的Degree值, Degree值越大,节点越大。选出排名前10位的基因作为Hub基因,分别是C-X-C趋化因子配体8(CXCL8)、Toll样受体2(TLR2)、细胞间黏附分子1(ICAM1)、选择素L(SELL)、趋化因子受体4(CXCR4)、细胞毒性T淋巴细胞相关抗原(CTLA4)、细胞分化抗原69(CD69)、双糖链蛋白多糖(BGN)、C-X-C趋化因子配体13(CXCL13)、金属蛋白酶抑制剂1(TIMP1)。将这10个关键基因建立Hub基因网络(见图6C)。

2.7 建立miRNAs-Hub基因调控网络

miRNA在基因表达的调控中起着多种作用。基于NetworkAnalyst数据库,使用Cytoscape构建miRNA-Hub基因调控网络。Hub基因及其相应的调控miRNA见图7。其中, miRNA (hsa-miR-335-5p)有5个靶基因(CXCL8、TLR2、ICAM1、CXCR4以及CD69); miRNA(hsa-miR-146a-5p)有4个靶基因(CXCL8、TLR2、ICAM1、CXCR4); miRNA(hsa-miR-92a-3p)有3个靶基因(ICAM1、CTLA4、CD69); miRNA(hsa-miR-155-5p)有3个靶基因(CXCL8、ICAM1、CTLA4); miRNA(hsa-miR-26b-5p)有3个靶基因(CTLA4、CXCL13、TIMP1); 另外,CXCL8、ICAM1是3个miRNA的共同靶标(hsa-miR-4426、hsa-miR-4462b和hsa-miR-4647);ICAM1、CD69是2个miRNA的共同靶标(hsa-miR-32-5p和hsa-miR-92b-3p);CXCL8、ICAM1是2个miRNA的共同靶标(hsa-miR-98-5p和hsa-miR-93-5p)。

3 讨 论

UC是IBD中的一种重要类型,其主要的特征是结直肠黏膜连续受累和弥漫型炎症。据认为, UC的发生是易感个体在环境因素作用下产生的先天性和适应性免疫功能障碍,导致的慢性、非受限制的炎症过程。这种炎症过程受遗传因素影响,并以免疫介导的方式持续存在。UC患者的病理学特点主要表现为结肠黏膜连续弥漫性慢性活动性炎症[15]。因此,探究发炎黏膜中的潜在生物标志物,对于获取新的治疗靶点尤为重要,同时对于后续临床诊断及治疗UC具有重要意义。

本研究利用GEO数据库中的GSE66407和GSE65114基因表达数据集进行了分析。首先,通过应用WGCNA和差异分析得到的差异基因取交集,获得了123个交集基因。其中, 117个基因表达上调,9个基因表达下调。然后对这些交集基因进行了功能富集分析,并构建了一个蛋白互作网络。这个过程筛选出了10个关键基因,包括CXCL8、TLR2、ICAM1、SELL和CXCL4等,这些基因在UC的发生和发展中可能起到重要的作用。

GO和KEGG功能富集分析结果显示,核心基因与白细胞迁移、白细胞介素、趋化因子活性以及胶原蛋白密切相关。涉及的信号通路包括PI3K-AKT信号通路、NF-κB信号通路、Toll样受体信号传导通路以及IL-17信号通路等,这些通路与炎症相关信号的激活和促炎因子的释放密切相关,提示其在UC的发病和加重过程中起着重要的作用。白细胞迁移参与IBD的进展,因此阻断白细胞向肠道迁移是控制疾病和缓解症状的主要策略[16]。研究[17]表明, PI3K-AKT细胞信号通路的异常被认为与UC的发病有关,对UC诱导的细胞凋亡和炎症以及长期UC诱导的结肠癌有很大影响。GO和KEGG富集通路分析结果表明,本研究中鉴定的关键基因可能通过上述途径参与UC的进展。

炎性趋化因子CXCL8是最早和研究最深入的趋化因子之一,可由多种细胞类型释放, CXCL8-CXCR1/2轴在炎症反应及肿瘤生长中起多效作用,包括中性粒细胞和粒细胞募集和活化到炎症部位,有助于消除病原体,在各种炎症性疾病中起重要作用。与正常对照组相比,活动性UC结肠黏膜中CXCL8表达增加已得到广泛证实, CXCL8表达上调与炎症程度相关[18-19]。TLR2作为模式识别受体(PRR)家族中的一员,能够感知病原相关分子模式(PAMP)和损伤相关分子模式。在机体免疫微环境中, TLR2信号通路对巨噬细胞、NK细胞、T细胞、肥大细胞等细胞的表型和功能产生影响,进而影响机体对病原体的清除能力和抗肿瘤能力[20]。ICAM-1是一种细胞表面糖蛋白和黏附受体,在许多细胞类型中充当黏附分子和信号受体,以引发炎症反应,启动炎症和愈合的消退,并调节肿瘤细胞的存活和传播[21]。SELL是一种细胞表面黏附分子,属于黏附归巢受体家族。该基因产物在白细胞与内皮细胞之间的结合和随后的滚动过程中起着必要的作用,从而促进其迁移到次级淋巴器官和炎症部位,同时还可以限制中性粒细胞的激活[22]。CXCL4可以改变单核细胞分化的轨迹,诱导包括Ⅱ类主要组织相容性复合体转录激活因子(CIITA)在内的关键转录调节因子介导的新型促炎和促纤维化表型[23]。

10个关键基因中,已有9个被报道[24-25]。与UC密切相关的基因包括CXCL8、TLR2、ICAM1、CXCR4、CTLA4、CD69、CXCL13、TIMP。此外,相关研究[26]通过随机森林算法和Limma分析得出, BGN可能是UC加重的基因特征。SELL与UC相关的报道较少见。研究[27]发现,SELL是牙周炎与非酒精性脂肪性肝病的共同关键基因。SELL还被发现是卵巢癌三级淋巴结构的相关特征基因,可作为预测卵巢癌预后和指导免疫治疗的生物标志物[28]。因此,进一步研究上述基因在UC中的作用,可以作为认识UC关键调控基因的新角度。

miRNA能够调节人类基因的表达,在自身免疫性疾病的发病机制中起着至关重要的作用[29]。糖皮质激素(GC)是治疗中度至重度活动性UC最有效的方法,研究[30]发现, miR-32-5p在GC耐药组中显著下调,热休克蛋白90β家族成员1(HSP90B1)受到了miR-642a-5p和miR-150-5p 的共同调控。相关研究[31]发现并验证了长链非编码RNA(lncRNA)MEG3的升高能够抑制UC大鼠中miR-98-5p的上调,并且通过调控miR-98-5p可以促进白细胞介素-10(IL-10)的表达。CHEN Z P等[32]认为,miR-146a-5p通过上调环指蛋白8(RNF8)和抑制Notch1/mTORC1通路减少肠上皮屏障损伤。PATHAK S等[33]发现,炎症介质诱导UC患者肠道肌成纤维细胞(IMF)中miR-155的表达,通过下调细胞因子信号抑制蛋白1(SOCS1)的表达。WU C P等[34]发现, Toll样受体4(TLR4)是hsa-miR-375的直接下游,其水平由has-miR-375负介导。hsa-miR-146a-5p、has-miR-155-5p、hsa-miR-335-5p、hsa-miR-32-5p、hsa-miR-98-5p、hsa-miR-93-5p及其靶基因已有与UC发生相关的报道[35], 本研究所得出的其余miRNA包括hsa-mir-92a-3p、hsa-mir-26b-5p、hsa-mir-4426、hsa-mir-4462b、hsa-mir-4647、hsa-mir-92b-3p及其靶基因可能在UC中发挥着关键的作用。

本研究使用WGCNA筛选差异基因,并且结合GSEA基因集富集方法和建立基因-miRNA调控网络,以探索UC潜在的调控机制,识别新的潜在生物标志物和通路,具有一定的创新性。基因-miRNA调控网络在UC病理生理学中起着重要作用。本研究的不足之处在于仅通过生物信息学分析微阵列表达谱,而没有通过相关实验进行验证以及探索Hub基因和miRNA调节UC的详细机制,因此,还需扩大临床研究样本和进一步实验以阐明机制。