缺血性脑卒中外周血差异miRNA挖掘与生物信息学分析

王成坦，唐娟，毛福秀，杨玉莹，袁珊，高蕊*

(1 石河子大学医学院/新疆地方与民族高发病教育部重点实验室，新疆 石河子 832000；2 石河子大学医学院第一附属医院神经内科，新疆 石河子 832000)

中风是危害人类健康的一类顽疾，其预后取决于缺血再灌注损伤诱发神经细胞的受损情况。研究表明，卒中后神经元死亡数目与缺血时间成正比[1]，因此，缺血性中风的早期诊断对于该疾病的治疗与预后意义非凡。然而，目前在临床应用方面尚未有一种可以单独应用于脑卒中诊断的外周血生物学标志物[1]。MiRNA(micro RNA)在人类脑组织和中枢神经系统中富集[2]，且对中风的发生、发展与预后发挥着重要作用[3]。本课题组前期研究发现，脑缺血再灌注后，低水平的miRNA-204可以减轻对其靶基因Nrn1的抑制作用，从而有助于形成适合神经修复的微环境[4]。在中枢神经系统疾病的发生发展中脑源性miRNA可以释放入血，引起外周血中miRNA水平改变[5]。

SONODA T等[6]通过miRNA芯片技术分析发现：外周血中有10个miRNA的表达变化与预测的中风风险之间存在显著相关性，其中有7个(miR-1268a、miR-1268b、miR-4433b-3p、miR-6089、miR-6090、miR-6752-5p和miR-6803-5p)在缺血性脑中风患者与正常受试者之间的表达具有显著差异性。这表明外周血中脑源性miRNA的表达变化与卒中的病理生理进程密切相关，可能具有潜在的诊断与治疗价值。随着生物信息学技术的发展与进步，利用生物信息学技术对前期脑卒中患者及动物卒中模型外周血中海量的实验性测序数据进行系统性的分析与挖掘，可能为寻找中风诊断标志物提供新的见解和方向。

大脑中动脉梗塞(MCAO，Middle Cerebral Artery Occlusion)模型是一种可以模拟脑缺血以及后续再灌注过程的动物模型，在研究缺血性中风中应用广泛。我们通过GEO(Gene Expression Omnibus)数据库下载了MCAO模型组与假手术组大鼠的外周血miRNA表达谱，应用R语言的limma包对测序数据进行miRNA差异分析，通过生物信息学方法预测其靶基因，并根据靶基因功能预测以及互作分析的连接度初步鉴定出核心靶基因。将上述核心靶基因在GEO数据库中的缺血性卒中患者外周血mRNA(message RNA)表达谱数据中进行鉴定，获得11个表达显著差异的关键基因。最后，结合miRNA-mRNA靶向数据以及关键基因的互作数据构建miRNA-mRNA调控网络，探讨上述调控通路在中风发生发展过程中可能发挥的生物学作用。

1 材料与方法

1.1 测序数据的获取与差异分析

在GEO数据库(http://www.ncbi.nlm.nih.gov/geo/)中，获取MCAO大鼠外周血miRNA表达数据集：GSE97532，其中包含3只MCAO大鼠与3只对照组大鼠。从GEO数据库中获取人类外周血mRNA表达谱数据集：GSE16561，涵盖39个缺血性中风病人的外周血样本与24个对照组受试者外周血样本，两组受试者的一般临床信息见表格一。

差异分析使用基于R语言编写的Limma软件包，R语言是一种基于S语言的开源程序语言，因开源开放且具有完整的开发者生态对其进行拓展和维护，因此被广泛应用于统计分析、绘图、数据挖掘等方面。使用Limma软件包进行差异分析的流程是：(1) 读取原始数据矩阵；(2) 矩阵清洗与标准化；(3) 分组与差异分析。差异筛选标准为：差异倍数≥1.5，P<0.05。

表1 缺血性卒中与对照组一般临床信息

1.2 差异miRNA的保守性分析与靶基因预测

应用miRbase数据库(http://www.mirbase.org/)鉴定差异miRNA在大鼠和人类序列之间的保守性，获取与人类序列高度同源的miRNA。应用miRDB软件(http://www.mirdb.org/)预测保守miRNA的靶基因，筛选标准为：Target Score>80。

1.3 核心靶基因筛选

将预测得到的miRNA靶基因上传至string数据库(https://string-db.org/)，预测靶基因之间的互作关系。将string数据库预测结果导入cytoscape软件(下载地址：https://cytoscape.org/)，分析并鉴定出连接度(degree)位于前20的核心靶基因。

1.4 关键基因的外周血表达水平鉴定

应用R语言的affy包对GSE16561芯片进行预处理，获取芯片表达量，表达量的归一化方法采用Z-Score法。提取20个核心基因在临床外周血样本中的表达信息，获得缺血性卒中患者外周血中表达水平出现显著改变的11个关键基因，差异分析使用 Student-T Test。

1.5 卒中后关键miRNA-mRNA表达网络构建

将差异miRNA与11个关键基因的靶向关系以及关键基因之间的相互作用关系合并后导入Cytoscape软件，构建关键miRNA-mRNA调控网络图。

1.6 关键基因的功能分析

使用R语言的ClusterProfile包对关键基因进行GO(Gene Ontology)与KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析。

2 结果

2.1 大鼠MCAO模型外周血差异miRNA的获取与保守性分析

我们应用R语言的limma包分析MCAO大鼠外周血miRNA表达谱GSE97532，得到15个差异miRNA，其中，上调miRNA有9个，下调有6个(图1A)。将差异基因按照表达量进行聚类，以鉴定差异基因的表达模式，直观展示差异miRNA表达水平(图1B)。应用miRbase数据库分析得到11个在人类和大鼠之间高度保守的miRNA：miR-346、miR-182、miR-194-5p、miR-139-5p、miR-128-2-5p、miR-24-1-5p、miR-191b、miR-107-5p、miR-383-5p、miR-328b-3p、miR-23b-5p，它们的序列在人和大鼠之间高度同源(图1C)。保守性高的miRNA往往提示其可能具有重要的生物学功能，因此在后续的分析过程中使用这11个miRNA。

图1 大鼠MCAO模型外周血差异miRNA获取与保守性分析

2.2 保守miRNA靶基因的预测与核心靶基因的筛选

MiRNA的生物学功能主要通过调控下游靶标基因的表达水平来实现。因此，我们采用miRDB软件预测miRNA的靶基因，设定严格的筛选标准，11条保守miRNA共预测到1467个靶基因。随后，应用string数据库，从基因共表达、基因融合、基因邻接、文本数据挖掘、实验证据证实等角度构建了差异miRNA靶基因之间的蛋白质-蛋白质互作网络(PPI network，Protein-protein Interaction Network)。在所有的1467个靶基因中，有1451个蛋白质被纳入本互作网络，组成含有1451个节点和7961对相关关系的复杂网络(图2)。

一般认为节点连接度(degree)越高的基因可能在互作网络中发挥更重要的作用，所以，我们应用cytoscape软件分析该互作网络并筛选出网络中连接度位于前20的核心基因(hub gene)。在筛选出核心靶基因之后，我们进一步的使用cytoscape软件从图2的PPI网络中提取出20个核心靶基因的互作关系网络(图3)。

图2 靶基因的相互作用网络图

颜色从黄色到红色代表着连接度(degree)由低到高，有互作关系的节点使用箭头相连，箭头所指为被调控关系图3 TOP 20 核心基因互作关系图

2.3 卒中患者核心靶基因外周血表达水平的鉴定

为了更进一步的明确核心靶基因与卒中疾病的关系，我们应用另一个mRNA的表达芯片GSE16561(包括39个中风病人外周血样本与24个对照组外周血样本)，检测20个核心靶基因的表达水平。结果发现有11个核心基因表达水平出显著差异(图4)。我们将这些在卒中患者外周血中显著变化且密切参与疾病发生发展互作网络的节点基因定为本研究的关键基因。绝大部分关键基因在卒中患者外周血中表达上调，仅有一个基因表达下调(图4)。这些在疾病中发挥重要生物学作用的关键基因，也具有成为卒中后外周血诊断标志物的潜能。

缺血性卒中组与对照组相比较，显著性检验使用Student-T Test，*P<0.05，**P<0.01，***P<0.001，****P<0.0001图4 卒中患者外周血关键基因的表达验证

2.4 卒中后miRNA-mRNA互作网络构建

在筛选并鉴定出差异miRNA的关键下游靶基因后，结合miRNA-mRNA靶向数据和关键基因蛋白质互作数据，应用cytoscape软件构建了卒中后miRNA-mRNA互作网络。从网络互作角度分析卒中后关键差异表达miRNA可能的作用。发现有6个miRNA(miR-182、miR-191-5p、miR-383-5p、miR-139、miR-128-2-5p、miR-194-5p)靶向调控11个关键基因，其中miR-182靶向调控了CTTN，CREB1，GRB2，FOXO3，FBXW7五个关键基因(图5)。

关键基因Notch1在互作网络中的节点连接度仅次于miR-182,Notch1以及与其具有相互作用的3个基因共同被筛选为本研究的关键基因(图5)。

红色V形为miRNA，绿色圆形为靶基因，线段代表相连的两节点之间具有互作关系图5 卒中后关键miRNA-mRNA互作网络

2.5 关键靶基因功能分析

我们应用GO和KEGG富集分析，对这11个关键基因的功能进行了预测和分析。GO分析结果显示：在生物学进程方面，造血调控、神经元死亡调控、对缺氧的反应等条目出现显著富集(图6A)；在分子功能方面，主要富集在与基因转录调控相关的转录因子激活与结合领域(图6B)，这些条目与卒中的发生发展密切相关。KEGG富集结果显示：PI3K-Akt信号通路、cAMP信号通路、Notch信号通路、MAPK信号通路等响应缺血性卒中的重要信号通路出现显著富集，且这几个信号通路之间存在相互作用，可以协同影响卒中的发生发展(图6C)。

图6 关键靶基因GO与KEGG富集分析结果

3 讨论

近年来缺血性卒中发病率不断上升，已成为造成中国成年人死亡的第一大主因[7]。不同于出血性卒中，缺血性卒中起病时没有典型的影像学特征表现，同时也缺乏可靠的分子诊断指标。在中风的相关研究中，有大量外周血表达变化的miRNA可能具有潜在的诊断和治疗价值，LONG G等人证实miR-30a，miR-126和let-7b可用作诊断缺血性卒中的可靠标志物[8]；MAKRIS K等发现miR-140-5p在卒中后外周血中表达显著上调，可能具备作为中风后抑郁的独立风险因素的价值[1]。LI等人报道，卒中后缺氧状态下上调的HIF-1a可以激活miR-182，抑制氧依赖性HIF-1a的降解酶PHD2和FIH的激活，由此形成正反馈促进HIF-1a的不可逆激活[9]。但是在临床应用方面，尚未有一种miRNA可以独立作为脑卒中的生物学诊断标志物[1]。

本研究挖掘得到的15个差异miRNA，其中有11个是大在大鼠与人类中高度保守的miRNA，而原始文献仅挖掘到7个差异miRNA，其中有5个与本研究相同[10]。原始文献应用的筛选软件是Affymetrix TAC，该软件属于芯片公司开发的软件，应用人数较少，分析较粗糙。本研究使用基于R语言的Limma包，它有以下几个优点：(1) 采用线性模型，适用范围广，在复杂的实验设计下也可以高效的处理数据；(2) 利用基因组数据的高度平行性，允许基因之间和样本之间存在不同程度的变异性，即使在小样本量的差异分析中也能保证分析的准确性；(3) 对表达谱数据整体进行分析处理，分析结果可以用于高层次基因表达特征分析、加权基因共表达网络分析等后续分析[11]。基于以上原因，limma包已成为应用最广泛的芯片分析方法[11]。本研究挖掘出来的miR-182、miR-128等关键miRNA，曾经被多项研究证实在缺血性卒中病人的外周血中出现显著表达变化[12-14]，可能在中风的诊断、治疗、预后等方面具有重要价值。

应用miRDB软件，11条高度保守的miRNA通过严格的筛选阈值，共筛选得到1467个靶基因。通过构建靶基因之间的PPI网络，分析网络连接度，我们找出了连接度位于前20的核心基因。通过分析缺血性卒中患者外周血mRNA表达谱数据，鉴定出卒中患者外周血中表达显著改变的11个关键基因。GO和KEGG富集分析结果显示，这些关键基因及其相关信号通路与缺血性卒中的发生发展密切相关。cAMP通路是关键基因高度富集通路之一，11个关键基因中的FOS、 CREB1、RHOA、CREBBP等基因都与cAMP通路有关。cAMP通路可以调节卒中后的神经再生、血管生成、神经炎症等生物学过程，是卒中后的重要调控通路[15]。CREB1是cAMP信号通路的关键蛋白，在癌症合并缺血性中风病人外周血中，CREB1基因显著高表达[16]。此外，Notch1在卒中后的脑血管内皮细胞中被激活，可以直接调节中风后的血管生成以及血脑屏障的完整性[17]。FOS基因在卒中梗死侧的周围持续高表达，并且在缺血性卒中病人的外周血中大量表达，这使其具备成为缺血性中风候选诊断标志物的潜能[18]。

结合miRNA-mRNA靶向数据、核心靶基因筛选、核心靶基因相互作用数据，本研究构建了卒中后的miRNA-mRNA表达调控网络，包含6个关键miRNA与11个核心基因。在这6个关键miRNA中，3个下调的miRNA:miR-194-5p、miR-182、miR-139调控了PPP2CA、FOS、Notch1、GRB2、CTTN、FOXO3、CREB1、FBXW7、ITPKB 9个功能基因的转录本，其中仅有ITPKB在卒中后表达下调，其余mRNA及其互补结合的miRNA的表达水平均呈反比关系。这与miRNA的作用方式有关，miRNA可以通过结合在靶基因的3’UTR区域沉默mRNA的表达[19]。3个上调的miRNA：miR-128-2-5p、miR-191-5p、miR-383-5p可能参与RHOA、CREBBP、FOS靶基因的调控。正如KEGG富集结果(图6C)所示，靶基因所参与的是一个复杂的调控网络，最终所体现的表达水平变化是各因素综合作用的结果。

结合靶基因所在的信号通路以及靶基因的表达谱数据，我们发现CREB1和Notch1在中风患者的血液样本中显著高表达(图4)，其相关通路——cAMP信号通路与Notch信号通路也被显著富集并具有重要的互作调控关系(图6C)。这提示我们，CREB1和Notch1以及靶向它的miR-182和miR-139可能在缺血性卒中的病理生理进程中发挥关键作用。前期研究表明，miR-182可作为中风早期诊断的外周血标志物[12]。miR-182可以通过调节靶基因的表达，影响中风后的组织损伤与修复过程。YI等[20]研究证实miR-182可以通过靶向调节ASPP家族的抑制因子来加重脑缺血再灌注损伤。而本研究发现，miR-182可能通过调控CREB1、GRB2、CTTN、FOXO3、FBXW7这5个关键基因影响cAMP信号通路、Notch信号通路、MAPK信号通路进而影响卒中后的病理生理进程。研究表明，miR-139可以通过直接靶向c-Jun调控缺血再灌注损伤后的炎症反应[21]。miR-139也可以通过靶向抑制Notch1进而抑制脑胶质瘤细胞的迁移侵袭能力和上皮-间质转化[22]。但是，miR-139靶向Notch1基因在脑卒中发生发展的作用尚未明确。我们证实了卒中后的关键作用因子miR-139和Notch1在外周血中的表达呈负相关，且miR-139和miR-182的靶基因作用途径高度相关，miR-182的靶基因GRB2可与miR-139的靶基因Notch1的调节剂Deltex直接作用，抑制MAPK信号通路的激活[23]。这些结果表明，miR-182和miR-139可能通过调控下游一系列靶基因发挥协同效应进而影响信号通路的互作，参与缺血性中风的病理生理进程，而其联合应用及作用机制有待更深入的探究。

综上所述，本研究通过生物信息学方法，鉴定了中风后差异表达的15个miRNA，通过同源性分析鉴定出11个人和大鼠之间高度保守的miRNA，预测并筛选了其靶基因。应用靶基因功能富集分析和靶基因之间的蛋白质互作网络解析，并与卒中患者外周血mRNA表达谱数据相结合，交集筛选出可能在中风的病理生理进程中起到关键作用的6个关键miRNA及其靶向的11个关键靶基因，构建了miRNA-mRNA表达调控网络，对该调控网络进行深度的功能分析。最终筛选出在缺血性卒中发挥关键作用的两对miRNA-mRNA靶向关系对：miR-182-CREB1和miR-139-Notch1，它们在缺血性卒中疾病中的联合作用机制和联合诊断价值值得更进一步的研究与关注。