【论著】基于基因表达综合数据库的缺血性卒中缺氧相关差异基因表达分析

摘要：目的 基于基因表达综合(GEO)数据库,采用生物信息学方法分析缺血性卒中的基因表达情况,结合缺氧相关基因,解析缺氧相关差异基因(HRDEGs)表达特征,筛选出关键基因,为深入了解缺血性卒中提供重要支撑｡ 方法 从GEO数据库下载GSE16561 和GSE58294 两个数据集,采用Python软件进行数据合并,利用Combat方法消除批次效应,同时保留疾病分组特征｡对消除批次效应前后的数据采用主成分分析法进行降维处理,并对缺血性卒中组和正常对照组人群进行组内相关系数(ICC)检测｡对合并及批次效应消除后数据集进行基因集富集分析(GSEA)及单样本GSEA,以名义P值(NOM P-val)< 0. 05且假阳性率P值(FDR P-val)< 0. 25作为标准,筛选出具有显著差异的基因集｡利用R软件对GSE16561和GSE58294数据集合并及批次效应消除后的缺血性卒中患者和正常对照者之间的差异表达基因进行鉴定[以log2 基因表达差异倍数(FC)的绝对值≥0. 58和矫正P值(Padj)< 0. 05作为筛选标准],并与在美国国立生物技术信息中心(NCBI)中获取的缺氧相关基因进行交集,得到HRDEGs｡对HRDEGs进行基因本体论(GO)和京都基因与基因组百科全书(KEGG)富集分析,使用STRING数据库构建差异表达基因的蛋白质互相作用网络,利用Cytoscape3 . 8 软件筛选排名前10 位的关键基因｡ 结果 ICC分析结果显示,消除批次效应后的GSE16561和GSE58294数据集中缺血性卒中和正常对照者的ICC分别为0. 94和0. 98,一致性极好｡GSEA结果显示,对GSE16561 和GSE58294 合并及批次效应消除后的新数据集中,34个基因集在缺血性卒中样本中显著富集,筛选出表达差异基因404个(均Padj < 0. 05),其中高表达基因354个,低表达基因50个｡与缺氧相关基因进行交集,获得64个HRDEGs｡GO富集分析表明,HRDEGs主要在囊泡腔､细胞质囊泡腔､分泌颗粒腔和特殊颗粒中显著富集,具有酰胺结合､肽结合､磷脂结合和酶抑制活性等分子功能,主要参与了细胞因子产生的正向调节､炎性反应的调节､对细菌来源分子的反应和对脂多糖的反应等生物学过程｡KEGG富集分析表明,HRDEGs主要在血脂与动脉粥样硬化､沙门氏菌感染､中性粒细胞胞外陷阱形成､核苷酸结合寡聚化结构域样受体信号通路､癌症中的蛋白聚糖､结核病和坏死性凋亡等通路上存在富集｡基于蛋白质互相作用网络,最终筛选出了10个关键基因,分别为ARG1(精氨酸酶1)､CASP1(含半胱氨酸的天冬氨酸蛋白水解酶1)､IL-1R1(白细胞介素1受体1型)､ITGAM(整合素亚基αM)､MMP9(基质金属蛋白酶9)､PTGS2(前列腺素内过氧化物合酶2)､STAT3(信号传感器和转录激活剂3)､TLR2(Toll样受体2)､TLR4､TLR8｡ 结论 通过生物信息学挖掘,本研究获得了10个缺血性卒中与缺氧相关的关键基因,可能为后续研究工作及诊断治疗提供了潜在靶点｡

缺血性卒中作为一种常见的脑血管疾病,是由各种原因引起的脑内血液循环障碍而导致脑组织的正常功能和生理结构发生改变｡鉴于神经细胞对缺血缺氧环境耐受性不佳,目前治疗的重点是快速进行再灌注治疗,即静脉溶栓和手术取栓治疗｡近年来,随着卒中的预防､救治､护理体系不断完善,个性化的预防干预措施相继涌现,我国卒中发病率及伤残调整寿命年逐年下降｡但溶栓药物和血管再通手术有着众多的禁忌证及治疗时间窗,导致最终受益人群受限｡目前,神经保护药物在治疗中无法起到关键性的作用,迫切需要新的预防策略及治疗缺血性卒中的有效靶点｡

近年来,生物信息学技术被广泛应用于发现复杂疾病相关的基因及机制通路,研究者通过尝试从不同方法､不同角度探索缺血性卒中中的重要靶点｡缺血性卒中所涉及靶基因众多､机制复杂,关键基因及通路仍需进一步挖掘｡

本研究尝试从分子生物学角度,分析缺血性卒中患者中相关基因表达的变化｡我们将两个独立芯片数据进行有效整合,对比了缺血性卒中组及正常对照组血液中基因差异表达情况,并综合了缺氧相关基因,采用生物信息学的手段,分析了缺血性卒中潜在的机制通路,筛选出缺血性卒中发生､发展过程中的关键基因,以为后续研究工作及诊断治疗提供有利支持｡

1 材料与方法

1. 1 数据来源

通过美国国立生物技术信息中心(National Center for Biotechnology Information, NCBI)平台上的基因表达综合(gene expression omnibus,GEO)数据库(https:/ / www. ncbi. nlm. nih. gov/ geo/),以“STROKE”为关键词检索有关缺血性卒中的人源血液样本基因表达数据集,以总样本数不低于20､样本数据包含正常对照组为条件进行筛选,同时排除外周血单个核细胞样本的数据集,最终选取GSE16561和GSE58294两组人源缺血性卒中患者相关外周血样本的基因表达芯片｡GSE16561数据集包含缺血性卒中组39例(男17例,女22例)和正常对照组24名(男10名,女14名),年龄均>40岁;GSE58294数据集包含缺血性卒中组69例(男36例,女33例)和正常对照组23名(男12名,女11名),年龄均> 50岁｡

本研究的两组数据均来源于GEO公开数据库,在研究中遵守GEO数据库使用规则｡

1. 2 数据集合并处理与验证

利用Python软件(版本号:3. 8. 10)对两组基因芯片数据合并,利用Combat 软件包(版本号:0. 3. 3)中Combat方法消除批次效应(同时保留疾病分组特征),并利用lets-plot包(版本号:4. 0. 1)绘制批次效应消除前后的箱图｡限于基因数量众多,为准确评估批次效应消除效果,应用主成分分析法(scikit-learn软件包,版本号:1. 3. 1)对消除批次效应前后的数据降维(累积方差解释率> 80% )处理｡然后利用Python 的pingouin 包(版本号0. 5. 3)按照缺血性卒中组和正常对照组分别对降维后的数据进行组内相关系数(interclass correlation coefficient,ICC)检测｡由于合并数据集来源的两组数据芯片属于数据样本不同的两组独立研究,选择ICC类型为ICC(3,k)｡根据ICC的值分别对消除批次效应前后正常对照组及缺血性卒中组数据进行一致性评价｡按照ICC检测标准,ICC的值介于0 ~ 1,< 0. 50表示一致性较差,0. 50 ~ < 0. 75 表示一致性中等,0. 75 ~ 0. 90 表示一致性较好,> 0. 90表示一致性极好｡最后,对消除批次效应前后的数据集分别利用lets-plot包作散点图｡

1. 3 基因富集分析(gene set enrichment analysis, GSEA)及单样本GSEA

利用Python语言GSEApy 软件包(版本号:1. 0. 6)对合并及批次效应消除后的数据集进行GSEA,选择来自人类分子特征数据库中的Hallmark基因集作为对照基因分析基因表达谱,并使用名义P值(nominal P value,NOM P-val)< 0. 05 且假阳性率P值(false discovery rate P value,FDR P-val)<0. 25作为标准,筛选出具有显著差异的基因集,同时使用GSEApy包中单样本GSEA算法对每个样本相应通路进行评分并使用Python 语言seaborn软件包(版本号:0. 12. 2)绘制热图｡

1. 4 获取缺氧相关差异基因(hypoxia relateddifferentially expressed genes,HRDEGs)

使用R软件(版本号:4. 3. 1)的limma 软件包(版本号:3. 56. 2)对合并及批次效应消除后的数据集进行差异基因表达分析,使用t检验分析并对P值矫正,以log2 基因表达差异倍数(fold change,FC)的绝对值和矫正P值(adjust P-value,P adj),即| log2 FC | ≥ 0. 58 且P adj < 0. 05 作为筛选高表达差异基因的标准,选出缺血性卒中患者与正常对照者之间的差异表达基因,并利用R 软件绘制火山图｡在NCBI(https:/ / www. ncbi. nlm. nih. gov / gene /)中检索“HYPOXIA”获取缺氧相关基因,从中挑选出缺血性卒中差异表达基因,得到HRDEGs,并使用R语言VennDiagram软件包(版本号:1. 7. 3)绘制韦恩图,随后使用Python 语言seaborn 软件包(版本号:0.12. 2)绘制HRDEGs 在合并及批次效应消除后的数据集中的聚类热图｡

1. 5 基因本体论(gene ontology,GO)和京都基因与基因组百科全书(Kyoto encyclopedia of genes andgenomes, KEGG)富集分析

使用R 语言clusterProfiler(版本号:4. 8. 2)､org.Hs. eg. db(版本号:3. 17. 0)､ggplot2(版本号:3.4. 2)软件包对缺血性卒中相关基因进行GO 及KEGG 富集分析并实现GO和KEGG结果的可视化｡

1. 6 HRDEGs蛋白互相作用(protein-protein interaction,PPI)分析

将HRDEGs 输入到STRING 数据库(https:/ /cn. String-db. org /),设置物种为“Homo sapiens”,构建PPI 网络｡形成的数据文件导入Cytoscape3. 8软件,利用软件中cytoHubba 插件从HRDEGs 中筛选网络中的重要节点｡本研究选择插件中最大集团中心度算法(maximal clique centrality,MCC)分析PPI,并将MCC 拓扑算法评分最高的前10 个基因定义为关键基因｡

2 结果

2. 1 数据集合并和批次效应消除处理及分析验证

对GSE16561 和GSE58294 这两个芯片数据集合并处理,初步合并后的数据集各样本基因表达情况如图1a｡利用主成分分析法对合并前后的数据进行降维处理(总计16 477 个基因数据信息,降维至40 个主成分,累积方差解释率> 80%),选用分析中的前2 个主成分绘制散点图如1b｡并对降维后的数据进行ICC 一致性检测,结果显示,合并后的数据集中缺血性卒中和正常对照者的ICC 分别仅为0. 42 和0. 76,两个数据集间存在明显的批次效应｡在保留疾病分组特征的前提下,利用Python 语言Combat 软件包的Combat 方法消除批次效应,处理后的数据集各样本基因表达情况如图1c;采取相同条件对消除批次效应后的数据进行降维处理,选用分析中的前2 个主成分绘制散点图如1d｡同时对降维后的数据进行ICC 一致性检测,结果显示,消除批次效应后的数据集中缺血性卒中和正常对照者的ICC 分别为0. 94 和0. 98,一致性极好,结合箱图及主成分分析散点图结果,批次效应消除｡

2. 2 合并及批次效应消除后数据集的富集分析

以Hallmark基因集中共计50个基因集为对照,分别对合并及批次效应消除后的数据集进行GSEA分析,发现其中34 个基因集具有显著差异(NOM P-val < 0. 05､FDR P-val < 0. 25),见表1｡

图2a ~ 2c分别示与代谢关系密切的缺氧､糖酵解和氧化磷酸化3个基因集在缺血性卒中组中富集｡对每个样本的相应通路进行单样本GSEA算法评分并用热图呈现结果如图2d所示｡

2. 3 差异表达基因的鉴定

对合并及批次效应消除后的数据集进行基因差异表达分析,结果绘制火山图(图3)｡共筛选出404 个差异表达基因,其中,表达上调的差异表达基因354个,表达下调的差异表达基因50个｡

2. 4 HRDEGs的筛选

在NCBI中检索“HYPOXIA”获取缺氧相关基因1 908个,与缺血性卒中差异表达基因的交集情况绘制韦恩图(图4),交集得到64 个HRDEGs,分别为PDK3(pyruvate dehydrogenase kinase 3,丙酮酸脱氢酶激酶3)､VASP(vasodilator stimulated phosphoprotein,血管扩张刺激磷蛋白)､IL-1R1(interleukin 1 receptor type 1,白细胞介素1受体1型)､ADM(adrenomedullin,肾上腺髓质素)､IL-1R2､PTBP3(polypyrimidine tract binding protein 3,多嘧啶束结合蛋白3)､RRM2B(ribonucleotide reductase regulatory TP53 inducible subunit M2B,核糖核苷酸还原酶小亚基M2B)､PDK4(pyruvate dehydrogenase kinase 4,丙酮酸脱氢酶激酶4)､TIMP2(tissue inhibitors of metalloproteinase 2,组织金属蛋白酶抑制因子2)､CASP1(cysteinyl aspartate specific proteinase 1,含半胱氨酸的天冬氨酸蛋白水解酶1)､MAPK14(mitogen-activated protein kinase 14,丝裂原活化蛋白激酶14)､ARHGAP26(Rho GTPase activating protein 26,Rho GTP酶激活蛋白26 抗体)､TLR4(Toll like receptor 4,Toll 样受体4)､CAMP(cathelicidin antimicrobial peptide,天冬氨酸抗菌肽)､STAT3(signal transducer and activator of transcription 3,信号传感器和转录激活剂3)､S100A12(S100 calcium binding protein A12,S100 钙结合蛋白A12)､HDAC4(histone deacetylase 4,组蛋白去乙酰化酶4)､SLC2A3(solute carrier family 2 member 3,溶质载体家族2 成员3)､MAP3 K2 (mitogen-activated protein kinase kinase kinase 2,丝裂原活化蛋白激酶激酶激酶2)､TXN(thioredoxin,硫氧还蛋白)､TNFSF10 (tumor necrosis factor superfamily member 10,肿瘤坏死因子超家族成员10)､ARG1 (arginase 1,精氨酸酶1)､MME (membrane metalloendopeptidase, 膜金属内肽酶)､DAAM2 (dishevelled associated activator of morphogenesis 2,形态发生紊乱相关激活因子2)､ANG(angiogenin,血管生成素)､TLR2､PFKFB3(6-phosphofructo-2-kinase / fructose-2,6-biphosphatase 3,6-磷酸果糖-2-激酶/果糖-2,6-二磷酸酶3)､HK2(hexokinase 2,己糖激酶2)､PYGL (glycogen phosphorylase L,糖原磷酸化酶L)､SLPI(secretory leukocyte peptidase inhibitor,分泌性白细胞肽酶抑制剂)､VMP1(vacuole membrane protein 1,液泡膜蛋白1)､HPSE(heparanase,乙酰肝素酶)､CREBRF(CREB3regulatory factor, CREB3 调节因子)､ANXA1(annexin A1,膜联蛋白A1)､ITGAM(integrin subunit alpha M,整合素亚基αM)､RGS2 (regulator of G protein signaling 2,G蛋白信号传导的调节因子2)､ORM1(orosomucoid 1,α-酸性糖蛋白1)､CCR7(C-C motif chemokine receptor 7,C-C 基序趋化因子受体7)､JAK2 (Janus kinase 2,Janus 激酶2)､HNRNPL (heterogeneous nuclear ribonucleoprotein L, 核不均一核糖核蛋白L)､ANXA3､PPP3R1(protein phosphatase 3 regulatory subunit B,alpha,蛋白磷酸酶3调节因子亚基1)､THBS1(thrombospondin 1,血小板反应蛋白1)､ABCA1 (adenosine triphosphate binding cassette subfamily A member 1,三磷酸腺苷结合盒亚家族A成员1)､PTGS2 (prostaglandin-endoperoxide synthase 2,前列腺素内过氧化物合酶2)､TLR8､DHX32 [DEAH-box helicase 32 (putative),DEAH 盒解旋酶32(假定)]､BIRC2(baculoviral IAP repeat containing 2, 含杆状病毒IAP 重复蛋白2)､WSB1(WD repeat and SOCS box containing 1,含Wd 重复和SoCs 盒蛋白1)､ANTXR2 (ANTXR cell adhesion molecule 2,ANTXR 细胞黏附分子2)､MTF1(metal regulatory transcription factor 1,金属调节转录因子1)､BMX (BMX non-receptor tyrosine kinase,BMX非受体酪氨酸激酶)､CEBPD(CCAAT enhancer binding proteindelta,CCAAT增强子结合蛋白δ)､VCAN(versican,多功能蛋白聚糖)､MMP9(matrix metallopeptidase 9,基质金属蛋白酶9)､AGO2(argonaute RISC catalytic component 2,argonaute RISC 催化组分2)､CASP5､SORL1(sortilin related receptor 1,索替林相关受体1)､SRPK1(SRSF protein kinase 1,SRSF蛋白激酶1)､PIK3CG (phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit gamma,磷脂酰肌醇-4,5-二磷酸3-激酶催化亚基γ)､ENTPD1 (ectonucleoside triphosphate diphosphohydrolase 1,外核苷三磷酸二磷酸水解酶1)､LRG1(leucine rich alpha-2-glycoprotein 1,富含亮氨酸的α-2-糖蛋白1)､SLC11A1 (solute carrier family 11 member 1,溶质载体家族11 成员1)､TCN2(transcobalamin 2,转钴胺素2)｡根据该64 个基因在合并及批次效应消除后的数据集中的变化情况绘制热图(图5)｡

2. 5 HRDEGs 的GO 和KEGG 富集分析

对获得的64 个HRDEGs,进行GO 富集分析｡如图6a ~ 6c 所示,HRDEGs 主要在囊泡腔､细胞质囊泡腔､分泌颗粒腔和特殊颗粒中显著富集,具有酰胺结合､肽结合､磷脂结合和酶抑制活性等分子功能,主要参与了细胞因子产生的正向调节､炎性反应的调节､对细菌来源分子的反应和对脂多糖的反应等生物学过程｡

KEGG 富集分析发现,HRDEGs 主要在血脂与动脉粥样硬化､沙门氏菌感染､中性粒细胞胞外陷阱形成､核苷酸结合寡聚化结构域样受体信号通路､癌症中的蛋白聚糖､结核病和坏死性凋亡等通路上存在富集(图6d)｡

2. 6 HRDEGs 的PPI

基于STRING 数据库,构建HRDEGs 的PPI(图7)｡

应用Cytoscape 3. 8 软件的cytoHubba 插件筛选出关联程度前10 的基因作为关键基因,分别为ARG1､CASP1､IL-1R1､ITGAM､MMP9､PTGS2､STAT3､TLR2､TLR4､TLR8,基于其表达蛋白质的PPI情况绘制表格｡见表2｡

3 讨论

缺血性卒中是一种机制复杂且靶点众多的疾病,许多基因在疾病过程中发生表达变化从而引起结构及功能上的改变,因而系统挖掘其潜在的关键基因以及分析其表达变化情况,为缺血性卒中的研究以及诊治方案可提供重要线索｡

本研究以血液为样本,不仅是由于血液为最普遍､易获取的样本,同时对于卒中而言,血液样本往往是全身效应的一个反馈,是全身病理生理状态的集中呈现形式,其在缺血性卒中的诊断､治疗及转归中均扮演着重要角色｡同时考虑到脑组织的低灌注是缺血性卒中的核心,将缺氧与缺血性卒中相结合能够极大地降低样本选择所带来的相关偏倚｡

本研究利用生物信息学手段,通过有效合并GEO 数据库GSE16561 及GSE58294 中缺血性卒中患者及正常对照人群血液样本的基因表达数据集,分别使用GSEA 及单样本GSEA 对合并及批次效应消除后的数据集进行富集分析,明确了缺血性卒中组及正常对照组外周血样本各通路活性状态,筛选出了缺氧､糖酵解､氧化磷酸化等34 组在缺血性卒中患者中显著富集的基因组｡同时我们从合并及批次效应消除后的数据集中获取了404 个差异表达基因,并与检索到的缺氧相关基因取交集,得到了64 个HRDEGs｡在对HRDEGs 进行富集分析中,GO 富集分析提示免疫应答及炎性反应在缺血性卒中中起到重要作用,KEGG富集分析结果同样表明差异基因主要富集于脂质动脉粥样硬化､免疫炎性反应及凋亡相关通路｡我们还利用PPI分析了上述64 个交集基因,从中筛选出了ARG1､CASP1､IL-1R1､ITGAM､MMP9､PTGS2､STAT3､TLR2､TLR4､TLR8 共10 个关键基因｡我们观察到筛选出的10个关键基因均与炎性反应和免疫应答相关,广泛作用于神经元､脑血管内皮细胞及神经免疫细胞｡

其中,TLR和IL-1R作为胞膜受体,在感知来自体内外组织异常中起着核心作用｡TLR 主要识别微病原体相关分子模式,同时也能感知细胞的非生理性破坏及细胞内容物的释放,而IL-1R 家族成员是内源性警报介质的受体,主要用于感知膜结构的完整性或细胞的适应性｡本研究结果显示,关基因中的TLR2､TLR4､TLR8 及IL-1R1 在缺血性卒中患者外周血中表达水平较高｡多项研究报道,缺血性卒中发生后,TLR2､TLR4､TLR8被激活,启动炎性反应｡反之,敲除上述基因则能够下调神经炎性反应,并改善预后｡然而有研究提示,虽然缺血性卒中患者血液中TLR2､TLR4 及TLR8 上调,但其表达量的高低与缺血性卒中患者的预后无关｡另有研究结果表明,TLR2､TLR4的激活在感染相关动脉粥样硬化中发挥重要作用｡同样,有研究指出,卒中能够诱导脑内皮细胞和胆碱能神经元中表达的IL-1R1 上调,激活IL-1,作用于血-脑屏障并促进中性粒细胞浸润,最终加剧了局部炎性反应损伤｡

PTGS2及STAT3 均能受到炎性反应诱导上调｡PTGS2上调会诱导前列腺素生成,进而促进炎性反应水平｡而JAK2 / STAT3 作为一种定位于细胞质的信号传导途径,受到脑缺血-再灌注损伤诱导激活,在调节卒中后炎性反应､细胞凋亡､氧化应激及血管生成中起到一定作用｡但在不同研究中,由于实验条件､药物及剂量､给药时间窗等不同导致最终研究结论不一致｡因此,STAT3 在缺血性卒中的具体功能及机制仍需进一步探讨｡

在缺血性卒中炎性反应的下游,关键基因ARG1及MMP9的表达量同样表现上调趋势｡缺血性卒中发生后,ARG1 不仅可以诱导中性粒细胞释放,也能够上调中枢神经系统中活化的小胶质细胞和巨噬细胞内ARG1 的表达,起到免疫抑制､调节炎性反应的作用,以控制卒中发展､促进损伤修复｡MMP9在缺血性脑损伤中同样起到双重作用,既能够降解细胞外基质,在卒中的早期破坏血-脑屏障,又可以为血管生成和胶质生成创造条件,有助于卒中后组织恢复及重塑｡

除此之外,脑缺血损伤还能够诱导炎症小体激活CASP1,使受损神经元发生焦亡;同时有研究证明,CASP1 能够促进氧糖剥夺神经元发生凋亡｡而ITGAM的产物CD11b +在中性粒细胞和单核细胞的黏附及吞噬中起重要作用,有研究表明,ITGAM高表达与缺血性卒中的动脉粥样硬化斑块不稳定和不良预后有一定关联｡但卒中后ITGAM表达调控的分子机制尚不清楚,仍需继续研究｡

然而本研究仍存在一定的局限性｡首先,血液作为研究对象,可受到全身各种因素的影响,尽管于临床诊断有提示作用,仍难免发生偏倚;其次,在研究中我们聚焦缺血性卒中,着重关注缺血后的血液变化,但缺血性卒中作为一类病因多样且机制复杂的疾病,疾病各亚型以及各阶段间的差异可能并不能被忽视｡

4 结论

生物信息学方法有助于把握缺血性卒中的研究方向,挖掘缺血性卒中的潜在治疗机制｡本研究筛选得到10个关键基因,表明缺血性卒中的发生发展与炎性反应､免疫联系密切,与ITGAM､MMP9､PTGS2､STAT3 等差异性表达基因的表达变化密不可分,筛选出的关键基因作为潜在的生物标志物及治疗靶点,可能有助于开发包括早期诊断和精确治疗在内的疾病管理的策略｡