基于生物信息学途径筛选缺血性脑卒中关键基因及药物预测

于 瑜，王钟兴

（中山大学附属第一医院麻醉科，广东广州 510080）

脑卒中作为全球第二大死亡原因，严重危害人类健康。全球每年约有550 万人死于脑卒中，每4个成年人中就有1人发病［1］。发病后遗留的身体残疾和神经功能障碍给中风幸存者的家庭及社会造成极大负担。在脑卒中病例中，缺血性脑卒中（以下简称脑缺血）约占70%，是脑卒中的主要类型［1］。目前，临床上脑缺血的主要治疗方法包括脑血管再通和神经保护。遗憾的是，这些措施收益有限。为了找到更有效的治疗方式以及治疗靶点，许多科学家从遗传学角度出发，寻找与脑缺血密切相关的基因，希望为治疗提供更多选择。例如：有研究人员利用高通量组学技术发现CLDN20、GADD45G、RGS2、BAG5 和CTNND2 可以作为小鼠和人类血液样本中脑缺血的血液生物标记物，协助脑缺血的临床诊断［2］。MAP2K4（丝裂原活化蛋白激酶4）是MAPK 通路的重要枢纽，MAP2K4 基因多态性rs3826392 可能在脑缺血后的炎症反应中发挥调节作用［3］。这些从基因层面解读脑缺血病理生理过程的研究对治疗方案的研发有一定启迪作用。因此，本研究借助生物信息学手段，通过对公共数据库的分析，寻找对缺血性脑卒中有重要作用的基因及相关通路，希望找到潜在治疗药物，能为缺血性脑卒中的进一步研究和后续治疗提供线索。

1 材料与方法

1.1 数据收集

大脑中动脉栓塞（middle cerebral artery occlu⁃sion，MCAO）是研究缺血性脑卒中的常用动物模型，可精确控制缺血时长。本研究选择GEO 数据库（https：//www.ncbi.nlm.nih.gov/geo/）中 研 究MCAO 小鼠表达谱变化的芯片集GSE98319。该数据集包含MCAO 小鼠和假手术组小鼠大脑皮层的检测结果。每组各3个样本。

1.2 差异表达分析

用R 语言的ggplot 包绘制6 个样本总体基因表达箱线图，并对芯片集进行多维度分析（multidi⁃mensional scaling，MDS）。Limma 包用于差异表达分析，根据测序平台注释信息及Ensemble 数据库对基因类型的注释信息，将探针ID 转换为相应的基因名称（gene symbol）并删除LncRNA 的数据。选择P<0.05 且|log2FC|>0.6 作为差异表达基因（dif⁃ferential expressed genes，DEG）筛选标准。差异分析结果以火山图形式呈现。用R 语言绘制热图展示数据集的聚类情况。

1.3 蛋白质互作网络构建及核心基因筛选

利用STRING（http：//www.string-db.org/）在线数据库对DEG 进行分析，预测可能在脑缺血发病过程中发挥重要作用的基因所编码蛋白质之间的互作关系。选择中度置信度（作用关系综合得分0.4-0.7）作为显著标准。将STRING 分析结果导入Cytoscape 3.7.1.软件进行可视化，使用Cytohubba和MCODE 插件对PPIN 进行分析。通过Cytohubba 计算出每个节点的度（Degree）［4］，Degree得分前10的基因构成的网络可视作核心子网络。通过MCODE（molecular complex detection）从蛋白质互作网络（protein-protein interaction network，PPIN）中筛出得分最高的3个功能模块亦视作核心子网络［5］。为了进一步缩小核心基因范围，将MCODE 评分最高模块的基因及Degree 得分前10 的基因进行整合，他们与DEG 中|log2FC|排名前20基因的交集即为核心基因（hub genes）。

1.4 小鼠大脑中动脉梗塞模型建立及核心基因表达水平测定

1.4.1 大脑中动脉栓塞模型建立及实验分组 3只6～8 周SPF 级雄性C57B/6J 小鼠购自湖南斯莱克景达实验动物有限公司，小鼠体质量为18～25 g，许可证号：SCXK（湘）2019-0004。本实验经中山大学实验动物管理与使用委员会同意。3%异氟烷诱导后将小鼠仰卧固定于操作台并戴上麻醉面罩，术中以1.5%异氟烷维持麻醉。颈部备皮消毒后颈正中切口，小心暴露右侧颈动脉鞘并将右侧颈总、颈外及颈内动脉分离出来，在三角分叉处远端结扎颈外动脉，近端穿一根线打活结，用动脉夹夹闭颈总动脉及颈内动脉，在颈外动脉中间剪一小口，将线栓插入颈外动脉后松开颈内动脉夹，将线栓缓慢沿着颈内动脉插入，线栓上的黑色标记点到达三角分叉处（插入约1 cm 左右），有明显阻力时停止。固定线栓，松开颈总动脉夹，缝合伤口。将小鼠放回笼内，侧卧，插栓侧朝上，维持体温。缺血30 min后经原切口拔出线栓，恢复血供，再灌注12 h 后留取小鼠双侧大脑皮质。小鼠缺血侧皮层为MCAO 组，对侧为对照组（CTRL组）。

1.4.2 标本收集以及qRT-PCR 验证基因表达变化 按照奕杉生物组织RNA 提取试剂盒说明书提取小鼠双侧皮层RNA，用翊圣生物Hifair® Ⅲ1st Strand cDNA Synthesis SuperMix for qPCR 试剂盒将RNA 逆转录成cDNA 后按照Hieff® qPCR SYBR®Green Master Mix 配置出实时荧光定量PCR（quanti⁃tative real-time polymerase chain reactionq，qRTPCR）体系进行检测。以Actin 为内参，用2-ΔΔCT相对定量法分析待测基因的转录水平。所用引物序列见表1。

1.4.3 统计方法 采用SPSS 17.0 统计软件进行统计分析。RNA 相对表达量为连续性资料，以均数±标准差（±s）表示，组间比较采用t检验（双侧），P<0.05认为差异有统计学意义。

1.5 GSEA算法进行功能富集分析

本研究选择GSEA C2 数据集中的KEGG 基因集（c2.cp.kegg.v7.2.）和C5 数据集中的GO 基因集（c5.go.mf.v7.2.entrez、c5.go.bp.v7.2.entrez、c5.go.cc.v7.2.entrez）为预设基因集，利用R 语言对差异分析结果进行GSEA 分析。FDR（false discovery rate）<0.05 作为富集通路的筛选标准。富集结果图由R 语 言的clusterprofiler 包［6］和GSEA 函数 进行绘制。

1.6 缺血性脑卒中药物预测

Connectivity Map（CMap）数据库包含1 309 种小分子化合物作用于多种细胞系引起的基因表达谱变化数据，可用于比较药物诱导的基因图谱与基因表达之间的相似性，并得到一个从-100 到100 的连通性评分：该评分大于0，表示化合物引起和上传基因相似的变化；评分小于0，表示该化合物引起和上传基因相反的变化，即该化合物可能对疾病有治疗作用［7］。由于CMap 将上传基因数目限制在150 个，本文将|log2FC|前150 上调DEG和全部下调DEG（100 个）分别上传至Cmap 在线数据库（https：//clue.io/）进行药物预测。连通性评分<-80 的小分子化合物作为本次预测的结果。

2 结果

2.1 缺血性脑卒中差异表达基因的筛选

GSE98319 数据集的箱线图提示各样本基因表达水平均衡（图1A），具有可比性。MDS（图1B）分析结果提示MCAO 组和假手术组能进行有效区分。图1C 显示GSE98319 芯片集差异表达分析的结果；不同变化方向的基因表达基因用不同颜色进行展示。将这些差异基因中lncRNA 删去，根据筛选条件，筛出521 个差异mRNA，其中421 个上调、100 个下调。图1D 为该数据集的热图，该数据集样本聚类效果较好，可信度较高。

2.2 缺血性脑卒中蛋白互作网络的构建以及核心基因的筛选

通过STRING预测521个DEG所编码蛋白质之间的相互作用关系。将预测结果导入Cytoscape 软件构建出蛋白质互作网络（图2）。此外，用Cyto⁃scape 的插 件CytoHubba 和MCODE 筛选出PPIN 中的核心子网络：度（Degree）排名前10的基因构成的子网络（图3A）以及MCODE 评分最高的子网络（图3B）。为了进一步缩小核心基因的范围，Degree 前10 的基因及构成MCODE 得分最高子网络的基因与DEG中上调倍数最高的20个基因求交集（图4），筛选出Drd4和Hcar22 个核心基因，表2 为他们的差异分析结果。

表1 qRT-PCR使用引物序列Table 1 Primers used for qRT-PCR

图1 GSE98319芯片数据集箱线图、MDS图、差异表达基因火山图以及热图Fig.1 Plots of GSE98319 gene chips

2.3 实时荧光定量PCR 验证核心基因

MCAO 组与CTRL 组各3 个样本，经统计分析发现脑缺血后Hcar2基因的相对表达量为（6.01±0.54）与对照组（1.00±0.09）相比有所上升，差异有统计学意义（t=15.79，P=0.000 094）。Drd4在脑缺血再灌注损伤组的相对表达量为（4.96±0.29）与对照组（1.00±0.02）相比，差异有统计学意义（t=23.74，P=0.000 019；图4）。qRT-PCR 结果与前述芯片差异分析结果一致。

表2 芯片数据核心基因差异分析结果Table 2 Results of analysis of microarray data for differential expression of hub genes

图2 差异表达基因的蛋白质互作网络Fig.2 PPIN of DEGs

图3 核心子网络Fig.3 Sub-networks of PPIN

图4 Hcar2 and Drd4 qRT-pcr验证结果Fig.4 Hcar2 and Drd4 expression level by qRT-PCR validation

2.4 GSEA功能富集分析结果

根据筛选条件，本次分析共富集到GO 分析中的12 个细胞组分（cellular component，CC）、16 种分子功能（molecular function，MF）和158 个生物过程（biological process，BP）以及11条KEGG 通路。图5和图6 分别显示了每一类富集结果中P值最小和富集分数最高的10条富集通路。

KEGG 结果显示细胞因子-细胞因子受体的相互作用、Toll 样受体信号通路、抗原加工呈递、胞质DNA 传感通路、RigⅠ样受体信号、MAPK 信号通路、自噬、自然杀伤细胞活性、NOD 样受体信号通路等与炎症反应密切相关的通路被富集，提示脑缺血再灌注损伤过程中炎症反应发挥了重要作用。GO富集分析得到了类似结果：对病毒的防御反应、病毒生命周期的变化、骨髓白细胞介导的免疫反应、对γ 干扰素和Ⅰ型干扰素的反应、对外部刺激反应的正向调节和白细胞的增殖以及细胞因子产生的正向调控相关基因在脑缺血后被显著富集，且全部处于激活状态，可能在脑缺血再灌注损伤中发挥重要作用；图5C、6C 显示，整合子复合体、高尔基体顺面膜囊、血液微粒、细胞外基质、分泌颗粒膜、质膜的外侧面、富含M-纤维胶凝蛋白的颗粒体、构成细胞外基质的胶原蛋白、嗜天青颗粒和双链断裂部位可能是脑缺血再灌注损伤中重要的细胞组分。图5D、6D 提示，脑缺血后脑组织在Ⅰ型干扰素受体结合、细胞因子活性、气味结合、趋化因子的活性、G 蛋白偶联受体结合、生长因子结合、清道夫受体的活性、生长因子的活性、丝氨酸水解酶活性、细胞因子受体结合等生物学功能方面发生了显著变化。在上述所有富集到的通路中，除高尔基体顺面膜囊的聚集为抑制外，其余通路均为上调或激活。

图5 KEGG、GO功能富集分析结果图Fig.5 Functional enrichment analysis results by GSEA algorithm

图6 GSEA功能富集结果Fig.6 Functional enrichment analysis results by GSEA algorithm

2.5 CMap药物预测结果

从CMap 下载预测药物结果，基于对药物的连通性评分进行排序和筛选。连通性评分<⁃80 的药物共有4 个：艾司洛尔Esmolol、甲巯咪唑Methima⁃zole、吐根酚碱Cephaeline、水仙环素Narciclasine。此4 种小分子化合物有可能成为缺血性卒中的潜在治疗药物（表3）。图7 为4 种小分子化合物的化学结构示意图。

表3 CMap预测结果Table 3 Prediction Results from CMap

图7 4种预测药物的简化结构式Fig.7 Structure of 4 potential drugs

3 讨论

本研究基于生物信息学途径筛选了脑缺血关键基因并预测了潜在药物。对结果进行分析发现构成核心子网络的基因大多与炎症反应相关，他们包括Toll 样受体（toll like receptor，TLR）家族（Tlr2、Tlr4）的成员、趋化因子CC-chemokine ligand（CCL）及相应受体CC-chemokine receptor（CCR）家族（Ccl5、Ccl2、Cxcl16、Ccr7），白细胞分化抗原Cd86、Cd14和Cd68等。其中Tlr2和Tlr4均已被证实在大鼠脑缺血再灌注损伤后表达水平显著升高［8］，同时甘草酸处理MCAO 大鼠可下调缺血脑内TLR2、高迁移率族蛋白B1（high mobility group box 1，HMGB1）及基质金属蛋白酶9（matrix metallopro⁃teinase-9，MMP-9）从而发挥神经保护作用［9］。脑缺血后TLR4/NF-κB 信号通路激活可能引起血小板聚集和星形胶质细胞激活从而加重神经损伤［10］。小鼠脑缺血72小时后调节T细胞表面CCR5表达增多，有利于减轻血脑屏障的损伤［11］。Georgakis等［12］基于全基因组关联分析发现较高水平CCL2有更高的缺血性卒中遗传倾向。这些结果均提示：炎症反应的调控对于脑缺血再灌注损伤的治疗有着巨大的意义。CD68、CD86、CD14、TLR2和TLR4都与巨噬细胞的功能密切相关。脑卒中患者血液中CD14highCD16+单核细胞的数量上升，这类单核细胞显著高表达TLR2 蛋白，同时CD14highCD16+单核细胞的高水平与高死亡率和脑卒中早期临床恶化相关［13］。结合本文分析结果，我们推测单核-巨噬细胞以及炎症反应的调控可能在脑缺血损伤中扮演了重要角色，可能是一个强有力的潜在治疗靶点。此外，ISG15-/-的MCAO 小鼠死亡率高，梗死加剧，神经系统恢复恶化，提示ISG15 可能作为内源性神经保护适应因子能减轻脑卒中损伤［14］。

有文献报道本文筛选出的核心基因Hcar2在小鼠脑缺血24 h后显著升高，48 h开始下降［15］。但其在脑缺血损伤中的作用以及机制仍不明确，另一核心基因Drd4编码多巴胺受体D4（dopamine receptor 4，D4R），能被多巴胺、肾上腺素和去甲肾上腺素激活参与情绪控制、视网膜节细胞的昼夜节律的调节［16］。其在脑缺血的病理生理过程中的作用亦有待进一步探索。艾司洛尔是一种选择性β1肾上腺素能受体拮抗剂，对脑缺血有一定保护作用：缺血前给予艾司洛尔能有效减轻双侧颈动脉闭塞引起的神经功能障碍［17］。在沙鼠脑缺血模型中，甲巯咪唑可诱导甲状腺功能减退以降低脂质过氧化、增加超氧化物歧化酶（superoxide dismutase，SOD1），使海马CA1 区神经元的死亡得到缓解［18］。艾司洛尔和甲巯咪唑的脑保护作用目前已相对明确，但具体作用机制尚未阐明。吐根酚碱和水仙环素是否对脑缺血损伤有保护作用则未见报道。

本文进行通路富集选用的GSEA 算法不需要事先筛选出差异基因，直接将全部基因与基因集进行比较，能保留表达变化不大但功能重要的基因，分析结果更可靠。本研究虽然通过实时荧光定量pcr 对Hcar2及Drd4的表达水平进行了验证，但未检测筛选出的药物的作用，仍须后续实验进行验证。Hcar2及Drd4以及4种潜在药物目前均未在脑缺血中得到充分的研究，可能是尚未发掘的新靶点。