基于生物信息学对心肌缺血再灌注损伤关键基因的筛选及实验验证*

王建茹， 李兴渊， 谢世阳， 程彦玲， 郭红鑫， 朱明军， 于 瑞

（河南中医药大学第一附属医院心脏中心，河南 郑州 450000）

缺血再灌注损伤是指由于栓子等原因堵塞动脉，导致组织或器官血液供应不足，缺血一段时间后，恢复血流导致组织缺血性损伤进一步加重，器官功能进一步恶化的综合征［1-2］。缺血再灌注损伤是许多疾病的主要病理因素，可发生在人体多个组织器官，如心、脑、肾、肝、肠等，临床可表现为心肌冬眠、心力衰竭、脑功能障碍、胃肠功能障碍、肾小管坏死、肝衰竭和多器官功能障碍综合征等［1-2］。心肌缺血再灌注 损伤（myocardium ischemia-reperfusion injury，MIRI）是缺血再灌注损伤最为常见的类型之一，急性心肌梗死经再灌注治疗后，导致心肌继发性损伤，再灌注所造成的心肌损伤可占最终心肌梗死面积的50%，已成为最常见的临床问题［3-4］。MIRI 病理过程复杂，涉及多种生物学途径，其具体机制尚未完全阐释清楚，临床防治仍具有挑战性，需不断探索和优化防治措施及方案。因此，积极探索MIRI 病理过程中潜在的分子机制及关键基因，发掘有效延缓或逆转MIRI的方法和措施，具有十分重要的临床意义，一直是心血管领域研究的重点和热点。

近年来，随着高通量技术及生物信息学的发展，他们在医学领域尤其是基础医学研究中得到了广泛的应用，可以帮助研究者揭示疾病的发生发展机制，挖掘疾病诊断和预后相关的生物标志物及调控网络，开发疾病新的防治措施或方法等。现阶段，高通量技术联合生物信息学分析的方法已被用于MIRI的研究中，为MIRI 的研究提供了强大的技术和方法支持，可帮助研究人员更全面地解析MIRI。诸多研究利用生物信息学分析的方法对高通量技术所获取的MIRI复杂数据进行了整合分析，发现MIRI病理过程中的铁死亡相关基因、竞争性内源RNA（compet‐ing endogenous RNA， ceRNA）调控网络、可能的关键基因、非编码RNA、潜在的防治小分子药物等［5-9］。总体而言，高通量技术和生物信息学在深化理解MIRI 病理机制，促进潜在药物靶点研发和诊断标志物的挖掘等方面做出了重要贡献。

因此，本研究拟从公共数据库中获取大鼠MIRI相关数据集开展相关研究（图1）。一方面利用稳健排序整合（robust rank aggregation， RRA）方法筛选稳健差异表达基因（differentially expressed genes，DEGs）；另一方面，利用替代变量分析（surrogate vari‐able analysis， SVA）包将各数据集合并后再筛选合并DEGs，取交集获取共同DEGs。随之，构建共同DEGs的蛋白相互作用（protein-protein interaction， PPI）网络，筛选关键基因。然后，构建MIRI 大鼠模型，检测关键基因的mRNA 和蛋白表达情况；并对关键基因介导MIRI 的情况开展文献回顾研究。最后，为了进一步揭示关键基因介导MIRI可能的机制对其开展基因集富集分析（gene set enrichment analysis， GSEA）。本研究将有助于理解MIRI病理过程中潜在的干预靶点和作用机制，可为其防治提供一定的参考。

Figure 1. Flowchart of research process. Limma： linear models for microarray data； SVA： surrogate variable analy‐sis； RRA： robust rank aggregation； DEGs： differen‐tially expressed genes； PPI： protein-protein interac‐tion； GSEA： gene set enrichment analysis； MIRI：myocardium ischemia-reperfusion injury. FC： fold change.图1 研究流程图

材料和方法

1 数据的下载和预处理

从GEO 数据库（http：//www. ncbi. nlm. nih. gov/geo/）和ArrayExpress 数据库（https：//www.ebi.ac.uk/biostudies/arrayexpress）中下载与大鼠MIRI 相关的数据 集GSE122020、E-MEXP-2098 和E-GEOD-4105。E-GEOD-4105 数据集分别由缺血30 min，再灌注2 d和7 d 两部分数据组成；GSE122020 和E-MEXP-2098数据集仅提取对照组（大鼠心脏前降支仅穿线不结扎）和实验组（大鼠心脏前降支结扎一定时间后再去除结扎线）的数据（表1）。依据平台/芯片的注释文件利用实用提取和报告语言（practical extracting and report language， Perl）脚本，将探针矩阵转为基因表达矩阵；利用avereps 函数对基因对应多个探针取均值；利用normalize Between Arrays 函数对数据进行矫正；基因表达数值较大者，取以2为底的对数。

表1 数据库中大鼠MIRI数据集的信息Table 1. Information of rat MIRI dataset in the database

2 稳健DEGs的筛选

利用微阵列数据线性模型（linear models for mi‐croarray data， limma）包，以|log2（fold change， FC）|>1和P<0.05 为阈值对预处理后的数据集的表达数据进行差异分析，筛选DEGs，获取每个数据集中基因的log2FC 值，并保存为TXT 文件；然后，利用RRA 包对TXT 文件进行整合分析，按照|log2FC|>1 和P<0.05为条件筛选稳健DEGs，并利用pheatmap 包分别将上调和下调前10的基因进行热图可视化。

3 合并DEGs的筛选

将各数据集合并为一个数据集，再利用SVA 包去除合并数据集内的批次效应，保存为“merged_da‐taset.txt”文件；然后利用limma 包，再以|log2FC|>1 和P<0.05 为条件筛选合并DEGs，并用热图进行可视化。

4 共同DEGs的获取及关键基因的筛选

利用Venn diagram 包将稳健DEGs 和合并DEGs取交集获取共同DEGs。将获得的共同DEGs 上传至STRING 数据库（https：//string-db.org/）中，选择“Mul‐tiple proteins”模式，种属设为“Rattus norvegicus”，互作得分设置为最高置信度≥0.400，进行蛋白互作分析，并将结果保存为tsv 格式。然后，将数据导入Cy‐toscape 3.7.2 软件，构建PPI 网络；利用cytoHubba 插件中的最大邻域组件密度（density of maximum neigh‐borhood component， DMNC）算法来筛选关键基因（即依据DMNC算法计算出的评分，按降序排序前5的基因）［10］。从“merged_dataset.txt”文件中提取关键基因的表达量，利用pROC包绘制关键基因的受试者工作特征（receiver operating characteristic， ROC）曲线以评价其诊断效能。

5 关键基因的验证

5.1 构建MIRI 大鼠模型 18 只6～8 周龄SD 大鼠［许可证号：SCXK（豫）2019-0002］随机分为对照组和实验组。2%戊巴比妥钠（40 mg/kg）腹腔注射麻醉并仰卧位固定大鼠，连接心电图，术区常规备皮消毒，气管插管，连接小动物呼吸机；于左侧第4、5 肋间切口，打开胸腔，暴露心脏，探查确定左心耳，于肺动脉圆锥与左心耳交界距离左心耳根部约3 mm 处用6-0号手术线结扎冠状动脉左前降支，肉眼见结扎线以下心肌颜色变苍白或发绀，心电图ST 段持续弓背抬高，说明心肌缺血成功；去除结扎线，心肌恢复红润且ST段下降至少1/2，提示再灌注成功。实验组以缺血30 min 再灌注2 h 来构建MIRI 模型，对照组前降支仅穿线不结扎［11］。该动物实验由河南中医药大学第一附属医院实验动物福利伦理审查委员会批准（批准号：YFYDW2020004）。

5.2 心肌梗死面积检测 利用2，3，5-氯化三苯基四氮唑（2，3，5-triphenyl-tetrazolium chloride，TTC）染色方法对心脏组织进行染色，以评价心肌梗死面积［11］。取下心脏后，PBS 洗涤3 次，利用大鼠心脏切片模具垂直于心脏长轴将其切成厚度为1mm的至少5 个切片，将切片置于2% TTC 染色液中，37 ℃避光孵育15 min，4%多聚甲醛固定，过夜拍照，梗死区为灰白色，非梗死区为砖红色。

5.3 血清心肌损伤相关指标水平的测定 腹主动脉取血后，室温静置1 h；4 ℃，3 000 r/min 离心15 min，分离血清。按照乳酸脱氢酶（lactate dehydroge‐nase，LDH）和肌酸激酶同工酶MB（creatine kinase MB isoenzyme，CK-MB）试剂盒说明书的操作步骤检测其水平。

5.4 关键基因mRNA 表达水平的检测 取两组大鼠适量心肌组织，加入研磨珠和RNA 提取液，用冷冻研磨仪充分研磨，提取总RNA 并测定其浓度，按照反转录试剂盒说明书合成cDNA 后进行RT-qPCR法，应用2－△△CT法对基因进行相对定量分析。引物序列由武汉赛维尔生物科技有限公司合成（见表2）。PCR 条件：95 ℃预变性30 s；95 ℃变性15 s；60 ℃退火延伸30 s，共40个循环。

表2 RT-qPCR引物序列Table 2. The sequences of the primers for RT-qPCR

5.5 Western blot 检测关键基因的蛋白表达情况取两组大鼠适量心肌组织，裂解，匀浆，离心提取总蛋白，BCA 法测定蛋白浓度，取80 µg蛋白上样，SDSPAGE 分离，转膜，封闭，孵育Ⅰ抗［MYC 原癌基因bHLH 转录因子（MYC proto-oncogene bHLH tran‐scription factor， MYC）、前列腺素内过氧化物合酶2（prostaglandin-endoperoxide synthase 2， PTGS2）、尿激酶型纤溶酶原激活物受体（plasminogen activator urokinase receptor， PLAUR）、胱天蛋白酶3（caspase-3， CASP3）和血红素加氧酶1（heme oxygenase 1，HMOX1）抗体，均1∶1 000稀释；β-actin 抗体，1∶2 000稀释］，4 ℃过夜，TBST 洗膜5 min×3 次，室温孵育Ⅱ抗（1∶5 000 稀释），洗膜，ECL 化学发光剂曝光显影，Image Lab 软件分析目的条带与内参照条带的吸光度，进行相对定量分析。

5.6 统计学处理 运用SPSS 21.0软件对进行统计分析，数据以均数±标准差（mean±SD）表示。两组间比较采用t检验。以P<0.05为差异有统计学意义。

6 文献回顾分析

分别以“关键基因”、“心肌缺血再灌注损伤”、“myocardial ischemic reperfusion injury”或“myocardi‐al ischemia reperfusion injury”为检索词，在中国知网（检索日期从1915 年1 月～2023 年7 月）和Pubmed（检索日期从1988 年1 月～2023 年7 月）中检索相关文献；下载并阅读文献，总结分析关键基因参与MIRI病理过程的情况。

7 关键基因的GSEA分析

以“merged_dataset. txt”文件中关键基因表达量的中位数为依据，将文件中的基因重新分为关键基因的高、低表达组。利用GSEA4.1.0 软件，选择c2.cp. kegg. v2023.1. Hs. symbols. gmt 作为参考基因集进行GSEA，以标准化富集得分（normalized enrich‐ment score， NES）的绝对数>1、名义P值（nominalPvalue）<0.05 和错误发现率（false discovery rate，FDR）<0.25为阈值筛选阳性基因集。

结 果

1 稳健DEGs的筛选结果

各数据集内样本间存在明显的差异（见图2A-1～2A-4），经过矫正后样本分布表现出了较好的一致性（见图2B-1～2B-4）。差异分析结果显示，以|log2FC|>1 和P<0.05 为阈值，E-GEOD-4105-2d 共筛选出172个DEGs（见图2C-1）；E-GEOD-4105-7d 共筛选出200个DEGs（见图2C-2）；E-MEXP-2098 共筛选出140 个DEGs（见 图2C-3）；GSE122020 共筛选出1 496 个DEGs（见图2C-4）。RRA 方法对各数据集进行整合，以|log2FC|>1 和P<0.05 为筛选条件，共鉴定出143 个稳健DEGs，其中上调114个，下调29个；前10个上调和下调稳健DEGs的热图见图2D。

Figure 2. The screening results of robust DEGs. A-1～A-4：box plots of samples in the E-GEOD-4105-2d， E-GEOD-4105-7d， EMEXP-2098， and GSE122020 datasets before preprocessing；B-1～B-4：box plots of samples in the E-GEOD-4105-2d， EGEOD-4105-7d， E-MEXP-2098， and GSE122020 datasets after preprocessing；C-1～C-4：volcano plots of DEGs in the EGEOD-4105-2d， E-GEOD-4105-7d， E-MEXP-2098， and GSE122020 datasets；D：heatmaps of top 10 upregulated and downregulated robust DEGs.图2 稳健DEGs的筛选结果

2 合并DEGs的筛选结果

各数据集合并为一个数据集后批次间存在明显的差异（图3A），经SVA 包去除批次效应后各样本呈现较好的一致性分布（图3B）。差异分析结果显示，以|log2FC|>1 和P<0.05 为阈值，在“merged_dataset.txt”文件中，共筛选出48 个合并DEGs，其中上调43个，下调5个（图3C）。图3D为前20个上调和下调合并DEGs的热图。

Figure 3. The screening results of merged DEGs. A：box plots of samples in the merged dataset before preprocessing；B：box plots of samples in the merged dataset after preprocessing； C：volcano plot of merged DEGs in the merged dataset；D：heatmap of the top 20 upregulated and downregulated merged DEGs.图3 合并DEGs的筛选结果

3 共同DEGs的获取及关键基因的筛选

分别将稳健DEGs 和合并DEGs 的上调和下调DEGs 取交集后，获得48 个共同DEGs，其中上调43个，下调5 个（见图4A、4B）。图4C 为共同DEGs 的PPI 网络，该网络包括221 条边和39 个上调节点。DMNC 算法鉴定出的关键基因包括MYC、PTGS2、HMOX1、CASP3）和PLAUR。为了进一步评价关键基因诊断MIRI 的效能，本研究基于“merged_dataset.txt”文件中的数据，绘制了关键基因的ROC 曲线，结果显示所有关键基因的曲线下的面积（area under the curve， AUC）值均大于0.8，除了PTGS2外其余4个关键基因的AUC 值均接近1，提示所有关键基因表现出了良好的诊断MIRI的性能（图4D）。

Figure 4. Acquisition of common DEGs and selection of hub genes. A： venn diagram of upregulated robust DEGs and merged DEGs；B： venn diagram of downregulated robust DEGs and merged DEGs； C： PPI network diagram of common DEGs； D： ROC curve plot of hub gene.图4 共同DEGs的获取及关键基因的筛选

4 关键基因的验证

本研究利用结扎SD 大鼠冠状动脉前降支缺血30 min 再灌注2 h 来构建MIRI 模型。对照组大鼠可见正常心电图，无ST 段、T 波、QRS 波的改变；实验组大鼠心电图可见ST段抬高，正常QRS波和T波消失，提示前降支结扎成功（图5A）。对照组未见心肌梗死，实验组可见明显的心肌梗死情况，血清LDH 和CK-MB 水平明显升高，提示MIRI 模型构建成功（图5B～D）。如图5E、F所示，与对照组相比，实验组大鼠心肌组织中CASP3、PTGS2、PLAUR 和MYC 的蛋白和mRNA 表达水平显著升高（P<0.05），HMOX1的蛋白和mRNA表达水平无显著差异（P>0.05）。

Figure 5. Validation results of hub genes. A： electrocardiogram results during MIRI induction in two groups of animals； B： TTC staining results of rat hearts in two groups； C： results of LDH activity in serum of rats in each group； D： results of CK-MB in serum of rats in each group； E： protein expression electrophoresis image of hub genes； F： relative protein expression lev‐els of hub genes； G： relative mRNA expression levels of hub genes. Mean±SD. n=3. #P<0.05 vs control group.图5 关键基因的验证结果

5 关键基因的文献回顾分析

通过文献回顾性分析，本研究发现PTGS2 的高表达加剧了MIRI 的程度；文献对MYC、HMOX1 和CASP3 参与MIRI 的报道存在不同的表达趋势，但未有报道PLAUR参与MIRI病理过程的研究（表3）。

表3 关键基因文献回顾的结果Table 3. Results of literature review on hub genes

6 PLAUR单基因GSEA结果

为了进一步揭示PLAUR 介导MIRI 病理过程可能的机制，本研究对PLAUR 开展了GSEA。结果显示，以|NES| >1、P<0.05 和 FDR<0.25 为阈值，PLAUR 高表达组共富集到19 个阳性基因集，主要涉及NOD 样受体信号通路、P53 信号通路、Toll 样受体信号通路、JAK-STAT 信号通路、细胞凋亡、趋化因子信号通路和自然杀伤细胞介导的细胞毒性等（图6A）；PLAUR 低表达组中富集到了9 个阳性基因集，主要涉及胰岛素信号通路、丙酮酸代谢、脂肪酸代谢和氨基酸降解等（图6B）。

Figure 6. Multiple GSEA enrichment plots of PLAUR. A： multiple GSEA enrichment plot of PLAUR high-expression group； B： mul‐tiple GSEA enrichment plot of PLAUR low-expression group.图6 PLAUR的多GSEA富集图

讨 论

MIRI 病理过程复杂，具体机制尚不完全清楚。已有研究表明，MIRI 的发生和发展涉及到多种病理机制和分子过程，如氧化应激、细胞凋亡、线粒体功能障碍、焦亡、非编码RNA、自噬和炎症反应等［27-28］。目前，MIRI 的治疗目标是减轻MIRI，保护心肌细胞免受损伤，主要的措施包括药物治疗（如抗凝血药、抗栓药、扩血管药、抗氧化剂和抗炎药等，旨在减轻MIRI 的程度）和介入治疗（如支架植入和血栓抽吸等，尽快恢复冠脉的血流，减少心肌缺血时间）等［29］。临床上，MIRI 现象仍无法避免，缺乏有效的治疗措施。积极探索MIRI 病理过程中的关键基因，具有重要的研究价值和意义，将有助于深入了解MIRI 的分子机制、寻找新的诊断和治疗标志物，可为其防治提供新的靶点和策略。

本研究一方面利用RRA 法对GSE122020、EMEXP-2098 和E-GEOD-4105 数据集中的DEGs 进行了整合分析，共鉴定出143个稳健DEGs，另一方面利用SVA 包将各数据集合并为一个数据集，共筛选出48 个合并DEGs；两者取交集后共获得48 个共同DEGs。随之，本研究构建了共同DEGs 的PPI 网络，并利用DMNC 算法鉴定出了5 个关键基因（MYC、PTGS2、HMOX1、CASP3和PLAUR）。这5个关键基因的AUC 值均大于0.8，提示其表现出了良好的诊断MIRI 的效能。动物实验研究结果显示，在大鼠MIRI心肌组织中，CASP3、PTGS2、PLAUR 和MYC 的蛋白和mRNA 表达水平升高，HMOX1 则无差异。文献回顾性分析提示，在MIRI 的病理过程中，无PLAUR 的相关研究，PTGS2 的高表达加剧了MIRI 的程度；MYC、HMOX1 和CASP3 在MIRI 病理过程中的表达存在上调、下调、无差异的情况。

MYC 又称c-MYC 属于Myc 基因家族，是一种多效性转录因子，可参与细胞的增殖、凋亡、分化、代谢、基因组稳定等多种细胞过程，与肿瘤的发生发展密切相关［30］。早期研究显示，MYC 的表达在猪MIRI的病理过程中没有明显变化［14］。随后，有研究报道MYC 在大鼠MIRI 及H9C2 心肌细胞缺氧/复氧模型中表达上调［12-13］；在小鼠MIRI 及小鼠心肌细胞缺氧/复氧模型中表达下调［15］。MYC具有调控细胞增殖和凋亡的双重作用，当受增殖信号刺激时MYC 表现出促进增殖效应，当受凋亡信号刺激时则表现出促凋亡效应［12］。PTGS2 是多种炎性前列腺素合成过程中的关键限速酶，其在大量的细胞类型中分布，可被细菌、脂多糖，促炎细胞因子和生长因子等激活，参与了肾、脑、心脏、肝脏、脊髓和肠道等脏器的缺血再灌注损伤的病理过；PTGS2 抑制剂可在一定程度上缓解上述脏器的缺血再灌注损伤的程度［31］。CASP3 是一种重要的半胱氨酸蛋白酶，属于半胱氨酸依赖性蛋白酶家族，在多种生理和病理过程中发挥着重要作用［32］。CASP3 作为一个关键的凋亡调控分子在MIRI 中扮演着重要角色，其异常表达和活化不仅调控心肌细胞的凋亡，还与炎症反应密切相关［32］。HMOX1 是一种应激蛋白和代谢酶，主要负责血红素的代谢过程，其可通过将血红素分解为碳一氧化物、游离铁和胆红素，具有抗氧化、抗炎、调节细胞信号和调节免疫等多种生理功能，在心血管疾病、代谢性疾病等多种疾病中发挥着重要作用［33］。实验证据显示，在MIRI 的动物和细胞模型中，HMOX1 的表达呈现出了高表达、低表达和无变化的情况［13-15］。本研究的生信分析结果提示HMOX1 在MIRI 中高表达，但动物实验验证结果却为HMOX1 在对照组和实验组间表达无差异。需要指出的是，HMOX1 在MIRI 病理过程中高表达不排除是代偿性升高的可能性。目前，有关HMOX1 调控MIRI 的作用机制仍需进一步的研究来探索。

PLAUR基因位于人染色体19q13.3 位置，可编码尿激酶型纤溶酶原激活因子受体蛋白。PLAUR是一种细胞膜表面受体蛋白，对于尿激酶型纤溶酶原激活因子（urokinase-type plasminogen activator， uPA）有较高的亲和力，能够特异性结合并激活uPA，继而将胞外的纤维蛋白溶解为可溶性产物，促进纤溶过程［34］。PLAUR在多种疾病的病理过程中均发挥了重要的作用，是一个具有广泛生物学功能的分子。研究显示，PLAUR与心肌梗死易感性相关，是动脉粥样硬化的有效诊断标志物，并表现出了一定的治疗价值；同时，PLAUR 的高表达还参与了脑、肾、肝移植、肌肉皮瓣等缺血再灌注损伤的病理过程［35-40］。目前，尚未见其与MIRI的关系，本文研究结果提示PLAUR的高表达介导了MIRI。为了进一步阐释PLAUR 介导MIRI 病理过程可能的作用机制，本研究开展了PLAUR 的GSEA，结果显示，NOD 样受体信号通路、P53 信号通路、Toll 样受体信号通路、JAK-STAT 信号通路、细胞凋亡和趋化因子信号通路等信号通路的激活与PLAUR 的高表达有关；胰岛素信号通路、丙酮酸代谢、脂肪酸代谢和氨基酸降解等途径的激活与PLAUR 低表达有关。总之，本研究首次创新性的提供了PLAUR 介导MIRI 的证据，认为PLAUR 可能是一个潜在的诊断和防治MIRI的靶点。

本研究尚存在一些局限性：（1）本研究虽运用RRA法和SAV包对多个数据集进行了整合分析和合并，但由于公共数据库中有关大鼠MIRI 研究的样本偏少，导致纳入分析研究的样本量有限。（2）本研究纳入的样本数据是利用基因芯片技术检测的。基因芯片技术在基因识别和注释等方面存在局限性，会导致某些基因可能无法被有效地识别和注释。此外，生物信息学分析的一些方法也存在一定的不足，比如在筛选DEGs 时，本文采用的是业界公认的|log2FC|>1 的标准；然而，|log2FC|=0.9 或0.8 的基因可能在MIRI 中也发挥着重要作用。以上情况可能会致使其他一些潜在的介导MIRI 病理过程的基因被忽略。（3）本文虽对挖掘出的关键基因PLAUR进行了文献回顾分析，动物实验验证以及单基因GSEA，但未针对探索出的PLAUR 介导MIRI 病理过程可能的分子机制进行功能验证，加之种属的差异，研究结果仍需要在后续的临床和基础研究中进行佐证。

总之，本研究利用生物信息分析的方法共挖掘出5 个关键基因（MYC、PTGS2、HMOX1、CASP3和PLAUR），通过文献回顾分析和实验验证后共鉴定出1 个（PLAUR）潜在的关键基因；PLAUR 可通过调控NOD 样受体信号通路、P53 信号通路、Toll 样受体信号通路、细胞凋亡、趋化因子信号通路、胰岛素信号通路和脂肪酸代谢等途径介导MIRI，结果可为进一步深入探索MIRI 潜在的防治靶点和分子机制提供了一定的借鉴和思路。