基于GEO 数据库和生物信息学分析筛选小鼠心肌缺血再灌注损伤相关的潜在枢纽基因

王建茹，彭广操，朱明军

河南中医药大学第一附属医院心内科，郑州 450099

急性心肌梗死（acute myocardial infarction，AMI）是冠状动脉粥样硬化性心脏病的一种急危重表现，发病率和死亡率均较高。近年来，随着人们生活方式的巨大变化，AMI 病例数量急剧增加，AMI 已成为我国居民住院和死亡的主要原因［1-2］。早期快速恢复冠状动脉血流可减少心肌梗死的面积，随着溶栓、经皮冠状动脉介入治疗等再灌注策略的发展，AMI 的死亡率有所降低。但大量证据表明再灌注会引起心肌的继发性损伤，其所造成的心肌损伤可占最终心肌梗死面积的50%，这种现象叫作心肌缺血再灌注损伤 （myocardium ischemia-reperfusion injury，MIRI）；MIRI 已成为AMI 经介入或溶栓治疗后最常见的临床问题［3-4］。MIRI主要表现形式多样，可进一步加重心脏结构和功能的损伤，最终引起多种不良心血管结局［5-6］。MIRI 的发生发展过程涉及炎症反应、氧化应激、钙超载、冠状微循环障碍、免疫反应、细胞死亡、长链非编码RNAs 表达水平变化等，但其详细的机制仍不清楚［7-9］。因此，积极探索MIRI 的分子机制，挖掘关键基因，对于MIRI的防治至关重要。

近年来，利用生物信息学技术对芯片基因表达数据进行挖掘已广泛地用于分析与疾病相关的差异表达基因（differentially expressed genes，DEGs），寻找关键基因，筛选与疾病诊断、治疗和预后相关的生物标志物等方面。DEGs有可能为探索MIRI病理过程的发展及防治提供有价值的线索。现阶段，研究人员可以从基因表达数据库（gene expression omnibus，GEO）等渠道获取芯片数据。但值得注意的是，在筛选DEGs 的过程中由于单项独立研究有限的样本数、样本的异质性、技术平台的差异、数据异常值及分析处理方法不恰当等因素的存在，研究结果可能出现一定的偏倚。目前，基于GEO 中基因表达芯片鉴定DEGs存在一些不足之处；因此，寻找更加有效的计算方法，可避免偏倚结果的出现，得到更加稳健的结果。稳健排序整合（robust rank aggregation，RRA）方法是一种可使多个数据集之间的偏差和误差达到最小的标准方法，用于各种排序基因列表的比较，可有效避免上述问题［10］。RRA 基于每个基因是随机排序的假设，如果一个基因在所有实验中排名高，用RRA 计算所得P值越小，则其成为DEG 的可能性越大。RRA 已广泛应用于多种肿瘤和非肿瘤疾病的多数据集的综合分析，并表现出了良好的效果。目前，尚没有关于RRA整合多个数据集分析MIRI的报道。

综上所述，本研究拟从GEO 数据库中获取小鼠GSE61592、GSE83472 和GSE160516 数据集，利用limma 包筛选各数据集中的DEGs，然后采用RRA整合分析3 个数据集的DEGs 筛选稳健DEGs；一方面利用clusterProfiler 包对稳健DEGs 进行基因本体论（gene ontology，GO）和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析，另一方面利用String 数据库构建稳健DEGs 的蛋白质-蛋白质相互作用（proteinprotein interaction，PPI） 网络；利用MCODE 和cytoHubba 插件来筛选PPI 网络中的子模块和hub 基因，再利用clusterProfiler 包对最重要的子模块基因和hub 基因进行GO 和KEGG 分析；最后，采用反转录- 定量聚合酶链反应（reverse transcriptionquantitative polymerase chain reaction，RT-qPCR） 技术检测小鼠MIRI 模型中hub 基因的表达情况。本研究旨在进一步探索MIRI 病理过程中潜在的干预靶点和分子机制，为其治疗提供新的策略和思路。

1 材料与方法

1.1 芯片数据的来源

本研究流程如图1 所示。从GEO 数据库中下载与小鼠MIRI相关的基因芯片数据集GSE61592［平台ID 为GPL6887，假手术（sham）组3 例，MIRI 组3例］、GSE83472（平台ID 为GPL6885，sham 组4 例，MIRI 组4 例）和GSE160516（平台ID 为GPL23038，sham 组4 例，MIRI 组12 例），其中GSE83472 的样本量信息是根据本次研究的需求提取的样本数，3 个数据集检测样本均为心肌组织。

图1 研究流程图Fig 1 Flow chart of the study

1.2 数据预处理和DEGs的筛选

利用perl脚本对表达数据进行注释，将探针矩阵转化为基因矩阵；利用impute 包对缺失值进行填充；利用avereps 函数对基因对应多个探针取均值；利用normalizeBetweenArrays 函数对数据进行矫正；若基因表达数值较大，则需转换为以2 为底的对数。利用limma 包，以|log2差异倍数（fold change，FC）|＞1 和P＜0.05 为阈值对3 个数据集的表达数据进行差异分析，筛选DEGs。

1.3 RRA法整合分析芯片数据

利用RRA 包对获得的DEGs 整合分析，按照|log2FC|＞1 和P＜0.05 为条件筛选稳健DEGs；然后利用pheatmap 包分别将P值排序前20 的上调和下调基因进行可视化。

1.4 稳健DEGs的功能富集分析

为了探究稳健DEGs 可能的功能作用，利用clusterProfiler包对其进行GO和KEGG富集分析。

1.5 PPI网络的构建

将稳健DEGs 上传至STRING 数据库，相互作用得分设置为最高置信度≥0.400，进行PPI 分析。利用Cytoscape3.7.2软件对PPI网络进行可视化。

1.6 模块分析和hub基因的筛选

将PPI 信息导入Cytoscape3.7.2 软件中，利用MCODE 插件并以默认参数（度截止=2，节点评分截止=0.2，最大深度=100，k 评分=2）为标准对PPI 网络进行模块分析，筛选最重要的子模块。CytoHubba插件可以通过其包含的12种算法分别对PPI网络中的每个节点（基因）赋予一定的分数，然后再依据基因得分将排序靠前的基因作为hub 基因。本研究的hub基因由12 种算法排序前100 的基因的重叠基因组成。然后，再利用clusterProfiler包对最重要的子模块基因和hub基因进行GO和KEGG富集分析。

1.7 Hub基因的验证

1.7.1 MIRI小鼠模型构建 18 只6～8 周龄的C57BL/6 J雄性小鼠购于郑州市惠济区华兴实验动物养殖场，许可证号：SCXK-（豫）-2019-0002。实验动物饲养于独立通风笼盒动物房，自由饮食和饮水。本实验已经河南中医药大学第一附属医院实验动物福利伦理审查委员会批准（批准编号：YFYDW2020004）。小鼠被随机分为sham 组和MIRI 组，每组9 只。戊巴比妥钠（40 mg/kg）腹腔注射麻醉并固定小鼠后，手术区备皮，消毒，于左侧第4、5 肋间切口，打开胸腔，撕开心包，轻压胸廓挤出心脏；探查确定左心耳，于肺动脉圆锥与左心耳交界距离左心耳根部约3 mm 处结扎冠状动脉左前降支，缺血30 min 后去除结扎线，关闭切口，缝合皮肤，24 h 后处死小鼠（即再灌注24 h）。Sham 组小鼠操作同MIRI 组小鼠，但不结扎冠状动脉左前降支。

1.7.2 氯化三苯基四氮唑染色 再灌注24 h 后，迅速摘取心脏，用磷酸盐缓冲液洗涤，从结扎线附近垂直于心脏长轴进行切片，将心脏平均切成5 片，然后置入2% 氯化三苯基四氮唑（triphenyl tetrazolium chloride，TTC）染色液（北京索莱宝科技有限公司）中37 ℃水浴15 min，结束后固定并拍照。TTC 染色后正常心肌组织呈红色，梗死组织呈苍白色。

1.7.3 RT-qPCR 各组小鼠取适量左心室组织（其中MIRI 组为心肌缺血危险区组织），加入TRIZOL（Invitrogen 公司，美国）进行裂解，振荡混匀，冰上放置5 min；加入氯仿，振荡静置离心；将离心管内的上清液转移至另一干净的EP 管中，加入异丙醇，振荡静置离心；去上清液，加入75%乙醇洗涤RNA，离心后吸掉多余液体，自然控干RNA；用无RNase水溶解RNA，轻轻吹打助溶，4 ℃静置10 min；用核酸蛋白浓度微量测定仪（Thermo Fisher 公司，美国）测RNA 浓度。按照反转录试剂盒（TOYOBO 公司，日本）说明书合成cDNA；按照2×Universal SYBR Green Fast qPCR Mix试剂（ABclonal公司，中国）的说明书配制反应液，在TL-988 荧光定量PCR 仪（西安天隆科技有限公司）上进行RT-qPCR；应用2-ΔΔCT对基因进行相对定量分析。引物由苏州金唯智生物科技有限公司合成（序列见表1）。PCR 反应条件：95 ℃3 min预变性，95 ℃5 s变性，60 ℃34 s退火延伸，共45 个循环。采用SPSS 21.0 软件对数据进行处理分析，数据以±s表示，2组间比较采用t检验，P＜0.05为差异有统计学意义。

表1 RT-qPCR引物序列Tab 1 Primer sequences for RT-qPCR

Continued Tab

2 结果

2.1 GEO数据集预处理和DEGs的筛选

差异分析结果显示，GSE61592 筛选出1 684 个DEGs，其中上调的832个，下调的852个；GSE83472筛选出55 个DEGs，其中上调的51 个，下调的4 个；GSE160516筛选出1 061个DEGs，其中上调的733个，下调的328个。图2为3个数据集的火山图。

图2 各数据集的火山图Fig 2 Volcano plots for each data set

2.2 RRA法筛选稳健DEGs

RRA 方法对3 个数据集进行了整合分析后，共鉴定出294 个稳健DEGs，其中上调的209 个，下调的85 个。图3 为P值排序前20 个上调和下调稳健DEGs 的热图。

图3 前20个上调和下调稳健DEGs的热图Fig 3 Heat map of the top 20 up-and down-regulated robust DEGs

2.3 稳健DEGs的功能富集分析

利用clusterProfiler 包、org.Mm.eg.db 包将筛选出的294 个稳健DEGs 进行功能富集分析，如图4所示。GO 分析结果显示，依据P＜0.05 确定了949个GO 条目；生物学过程（biological process，BP）条目为867 条，主要涉及细胞的迁移、趋化作用、对外界刺激反应的正性调节、调节肿瘤坏死因子超家族细胞因子的产生、细胞因子介导的信号通路等；分子功能（molecular function，MF）条目为62 条，主要涉及受体配体活性、细胞因子受体结合、细胞因子活性、趋化因子活性、细胞因子受体活性、Toll 样受体结合等；细胞组分条目（cellular component，CC）为20 条，主要涉及细胞外基质、线粒体外膜、肌质网、溶酶体等。KEGG分析结果显示，依据P＜0.05 共映射出5 条信号通路，包括细胞因子受体相互作用、病毒蛋白与细胞因子及细胞因子受体的相互作用、阿米巴病、白介素-17 和核因子κB （nuclear factor kappa B，NF-κB）信号通路。

图4 稳健DEGs的功能富集分析Fig 4 Functional enrichment analysis of robust DEGs

2.4 稳健DEGs的PPI网络构建和模块分析

利用STRING 数据库对筛选得到的稳健DEGs 构建PPI网络（图5A）。该网络包括1 248 条边和239 个节点，其中上调的174个，下调的65个。模块分析共筛选出14 个子模块，其中模块1 最重要（图5B）。模块1 中的稳健DEGs 主要富集在急性炎症反应、炎性细胞的迁移、细胞趋化、细胞黏附分子结合、趋化因子受体活性、病毒蛋白与细胞因子及细胞因子受体的相互作用、NF-κB 信号通路、趋化因子信号通路等（图5C、D）。

图5 稳健DEGs的PPI网络构建和模块分析Fig 5 PPI network construction and modular analysis of robust DEGs

2.5 Hub基因的筛选

共获得17 个hub 基因（图6A、B），相关信息详见表2。功能富集分析的结果显示，17 个hub 基因主要富集在白细胞和中心粒细胞的迁移、淋巴细胞趋化性、对外界刺激反应的正性调节、趋化因子的活性、趋化因子受体结合、细胞因子受体结合等功能；涉及病毒蛋白与细胞因子及细胞因子受体的相互作用、细胞因子受体相互作用、趋化因子信号通路、NF-κB信号通路等通路（图6C、D）。

表2 Hub基因的相关信息Tab 2 Information about the hub genes

图6 Hub基因的筛选及其功能富集分析Fig 6 Screening of hub genes and analysis of their functional enrichment

2.6 Hub基因的验证

由图7A 可知，与sham 组相比，MIRI 组小鼠经过缺血再灌注处理后，心肌明显出现了梗死，说明MIRI 模型构建成功。如图7B 所示，与sham 组相比，MIRI 组小鼠心肌趋化因子（C-C 基序） 配体4［chemokine（C-C motif） ligand 4，Ccl4］、Ccl6、Ccl7、趋化因子（C-X-C 基序）受体4［chemokine（C-X-C motif） receptor 4，Cxcr4］、趋化因子（C-C基序） 受体2 ［chemokine（C-C motif）receptor 2，Ccr2］、信号调节蛋白β1（signal-regulatory protein β 1，Sirpb1）、低亲和力免疫球蛋白γ Fc 区受体Ⅱb（low affinity immunoglobulin gamma Fc region receptor Ⅱb，Fcgr2b）、白细胞表面抗原Cd53（leukocyte surface antigen CD53，Cd53）、花生四烯酸5-脂氧合酶激活蛋白（arachidonate 5-lipoxygenase activating protein，Alox5ap）、髓样分化初级反应基因88（myeloid differentiation primary response gene 88，Myd88）、巨噬细胞清道夫受体1 （macrophage scavenger receptor 1，Msr1）、 基质金属肽酶14（matrix metallopeptidase 14，Mmp14）、髓样细胞上表达的触发受体2 （triggering receptor expressed on myeloid cells 2，Trem2）、桩蛋白（leupaxin，Lpxn）的mRNA 表达量均上调，而低亲和力免疫球蛋白γ Fc 区受体Ⅲ（low affinity immunoglobulin gamma Fc region receptor Ⅲ，Fcgr3）、补体C1q亚组分亚单位B（complement C1q subcomponent subunit B，C1qb）、去整合素和含金属蛋白酶结构域蛋白（a disintegrin and metalloproteinase domain-containing protein 8，Adam8）的mRNA表达未见明显的差异。回顾文献发现，目前仅Fcgr3、C1qb、Trem2、Lpxn、Cd53、Alox5ap、Sirpb1、Fcgr2b在MIRI 与正常组织中的差异表达未被报道。

图7 Hub基因的验证Fig 7 Validation of the hub genes

3 结论

本研究结合微阵列基因表达谱和生物信息学方法，探索了MIRI 病理过程中可能的新的核心基因及其潜在的分子机制。研究首先利用limma 包对小鼠GSE61592、GSE83472 和GSE160516 数据集中的DEGs 进行了筛选。其次，利用RRA 法对DEGs 进行合并分析，共鉴定出209 个上调和85 个下调稳健DEGs。然后，对稳健DEGs的PPI网络进行功能模块分析，共发现14 个子模块，其中模块1 最重要；同时，鉴定出了17个hub基因，而且这些基因大部分都在模块1 中。最后，本文对稳健DEGs、模块1 中的基因和17个hub基因分别进行了功能富集分析，结果发现这些基因大都是通过调控炎性细胞的迁移、趋化因子和细胞因子及其受体活性、细胞因子受体相互作用、趋化因子信号通路、Toll 样受体信号通路、NFκB 信号通路等生物学功能和通路，来介导MIRI的病理过程。众所周知，MIRI 是一个复杂的事件，涉及炎症细胞（如中性粒细胞、单核/巨噬细胞、淋巴细胞） 和关键激活信号（如细胞因子、趋化因子）等［11］。研究报道，炎症信号通路Toll 样受体和NFκB 信号通路以及两者之间的交互作用可加剧MIRI的严重程度［12-14］。此外，MIRI可强化趋化因子的表达，这些趋化因子可诱导炎性细胞的活化和浸润，进而进一步产生强烈的炎症反应，而抑制趋化因子可促进损伤心肌的愈合［11，15］。以上相关文献的报道从侧面证实了本研究富集分析结果的科学性。

为了进一步确证生物信息学分析预测的结果，本研究构建了小鼠MIRI 模型并用RT-qPCR 方法检测了17 个hub 基因mRNA 的表达情况。结果显示，在MIRI 过程中Ccl7、Sirpb1、Fcgr2b、Cxcr4、Cd53、Alox5ap、Ccl6、Myd88、Msr1、Ccl4、Mmp14、Ccr2、Trem2、Lpxn的表达上调，而Fcgr3、C1qb、Adam8未见明显差异。既往其他研究也证实了在MIRI 的过程中Ccl7、Cxcr4、Ccl6、Myd88、Ccl4、Msr1、Ccr2、Mmp14的表达上调［16-22］，这与本研究的结果吻合，从侧面佐证了本研究结果的可靠性和科学性。Adam8参与了MIRI 的病理过程，但在再灌注0.5 h 和24 h 时其表达并没有增加，48 h 后才明显增加［23］。目前，Fcgr3、C1qb、Trem2、Lpxn、Cd53、Alox5ap、Sirpb1、Fcgr2b在MIRI中的差异表达情况，均未被其他研究报道；但Fcgr3和C1qb在本研究中的生物信息学分析预测结果与动物实验的验证结果不一致。综上所述，Trem2、Lpxn、Cd53、Alox5ap、Sirpb1、Fcgr2b这6 个hub 基因的高表达在MIRI 中发挥着重要作用，将有助于阐释MIRI可能的分子机制。

Alox5ap位于染色体13q12-13，与心肌梗死、动脉粥样硬化等心血管疾病的易感性密切相关，其编码的5-脂氧合酶激活蛋白，可将花生四烯酸催化为白三烯，进而参与炎症反应介导多种心血管疾病［24］。Fcgr2b作为唯一的抑制性Fc 受体，在调控炎症和免疫方面发挥着重要作用，与自身免疫性疾病的易感性密切相关［25］。Fcgr2b广泛表达于巨噬细胞、树突状细胞等免疫细胞，可抑制激活FcγRs诱导的吞噬前信号，进而抑制吞噬作用［25］。研究［26］显示，高效的吞噬作用可以及时清除MIRI过程中凋亡的心肌细胞，否则，将引起炎症反应、细胞进一步凋亡等后续损伤。因此，Fcgr2b可能通过调控吞噬作用来参与MIRI的病理过程。Trem2为免疫球蛋白超家族受体的成员，在小胶质细胞、巨噬细胞、树突状细胞、髓系细胞等细胞中表达，可与相对分子质量为12 000的跨膜免疫接头DNAX 激活蛋白（DNAX-activating protein of 12 kDa，DAP12）的酷氨酸活化基序结合，通过激活下游信号通路，在免疫炎症等反应中发挥作用［27-30］。Trem2已被报道通过调控细胞焦亡及凋亡、炎症反应等，参与了脑、肝、肾、肺缺血再灌注损伤的病理过程［27-30］。由于各器官缺血再灌注存在相似的病理生理过程，因此Trem2可能也介导了MIRI。研究［31-32］显示，B 细胞被募集到了MIRI 的心肌组织，其加剧MIRI 的机制可能与其分泌IgM 进而激活补体系统和炎性细胞有关；同时，被募集到MIRI 心肌组织的还有T 细胞尤其是CD4+T 细胞，其主要表现为Th1 促炎症表型，通过释放γ 干扰素等，加剧MIRI。Cd53为四次穿膜蛋白（tetraspanin）超家族的成员，高表达于T 细胞、B 细胞等免疫细胞表面，在免疫系统中发挥着重要作用，主要表现为调控免疫细胞的黏附、迁移、增殖、活化以及信号转导等［33］。Lpxn是桩蛋白（paxillin）家族的成员，表达于白细胞细胞系、B 细胞等细胞，可调控整合素的活性，在B 细胞免疫反应中发挥着不可或缺的作用［34］。Sirpb1是免疫球蛋白超家族的信号调节蛋白成员，胞内段缺乏信号转导基序，但可通过与DAP12结合，转导活化信号，参与细胞增殖、分化、免疫功能等［35］。值得注意的是，Sirpb1与B 细胞受体通路相关，参与了B 细胞的增殖［36］。也就是说，Cd53、Lpxn、Sirpb1可能通过调控免疫机制来介导MIRI 的病理过程。综上所述，我们尚未发现有关这6 个hub 基因在MIRI 中作用的相关报道，但回顾文献的结果提示这些hub 基因表现出了介导MIRI 病理过程的可能性，为探索它们在MIRI 中的作用机制提供了一定的依据。

本研究运用生物信息学的方法发现了在小鼠MIRI 过程中的一些hub 基因，但仍有一些局限性。虽运用RRA 法分析了多个数据集，但由于基因芯片数据的可获得性有限以及每个数据集中满足研究条件的样本较少，最终导致本研究纳入的样本量有限。生物信息学预测结果虽进行了实验验证，但由于种属的差异，结果仍需要进一步在临床实践过程中进行佐证。总之，本研究利用基因表达谱和生物信息学相结合的方法共挖掘出17个hub基因，通过回顾既往文献和实验验证后共鉴定出6个潜在的hub基因（Fcgr2b、Alox5ap、Trem2、Sirpb1、Cd53、Lpxn），研究结果可为进一步探讨MIRI 的分子机制和治疗靶点提供新的思路和切入点。