基于生物信息学方法筛选缺血性心肌病生物标志物的研究

张开健，胡康，张步春

缺血性心肌病(ischemic cardiomyopathy，ICM)为冠心病晚期阶段和特殊类型，主要发病机制与冠状动脉粥样硬化有关，引起心肌变性、坏死和纤维化，随着疾病进展，容易发展为心力衰竭。而ICM心力衰竭具有起病急、病情进展快、致死率高等特点，随着人口老龄化，ICM发病率呈升高趋势[1-2]。生物信息学是研究各种组学数据，包括转录组学、蛋白质组学等生物学数据的一门新兴交叉学科，研究方法主要包括对生物学数据的搜集(收集和筛选)、处理(编辑、整理、管理和显示) 、利用(计算和模拟)及分析[3]，近年来实践表明生物信息学技术在筛选疾病生物标志物中有极大的应用价值，其对疾病的诊断、治疗以及预后有极大的意义，使人们对疾病的认识更加全面且深刻。目前许多预测心血管疾病生物标志物已经被评估，且成功地应用于心血管疾病的风险预测[4]，但在心肌缺血方面还需要进一步研究确定新的生物标志物。本研究基于生物信息学方法，通过分析ICM患者心肌组织的高通量测序芯片，初步筛选心肌缺血相关的潜在生物标志物，以期为ICM诊疗提供依据。

1 材料与方法

1.1 数据来源 研究数据来自美国国家生物信息中心NCBI(www.ncbi.nlm.nih.gov)中的GEO数据库，利用R语言(R4.1.2)直接从GEO数据库中下载数据集GSE26887，该数据集用取自患者心脏左心室活检组织，基于GPL6244检测平台，用Affymetrix GeneChips Human Gene 1.0 ST芯片检测基因表达，包含5例健康对照、7例缺血后扩张型心肌病合并糖尿病、12例缺血后扩张型心肌病不合并糖尿病患者的基因的RNA表达谱。本研究设置5例健康样本为对照组(GSM662158—GSM662162)、12例ICM患者为实验组(GSM662179—GSM662181)，2组患者的年龄、性别、吸烟、高血压、体重指数(BMI)等一般资料比较差异无统计学意义，且排除了重要影响因素糖尿病对研究的影响。

1.2 数据的获取与前期处理 基于R语言(R4.1.2)的GEOquery包获取GSE26887数据集，通过GPL6244平台对应的hugene10sttranscriptcluster.db包对数据集基因芯片探针进行基因名称注释，获取基因名称、基因表达值，去除探针未对应基因名的数据及一个基因名对应多个探针的数据(仅保留最大表达量探针对应基因名)，去除7例ICM合并糖尿病数据集。

1.3 DEGs的获取与可视化 基于R语言limma包对上述数据进行差异表达分析，获得DEGs及LogFC、t值、P值、adj.P等，由于存在DEGs基数较大，部分DEGs差异倍数(fold change，FC)较小、差异无统计学意义等因素，设定筛选条件|LogFC|>1且adj.P<0.01进一步筛选DEGs，并基于R语言ggplot2包、pheatmap包分别对其绘制火山图、热图进行可视化处理。

1.4 DEGs的GO和KEGG富集分析 GO(gene ontology，GO)富集分析可以从基因的生物学过程(biological process，BP)、细胞组分(cellular component，CC)及分子功能(molecular function，MF)等方面进行基因注释，KEGG(kyoto encyclopedia of genes and genomes，KEGG)富集分析可以对DEGs的信号通路方面进行富集分析，以挖掘疾病相关基因的生物学通路，基于R语言clusterProfiler、org.Hs.eg.db包行DEGs的GO和KEGG富集分析，使用ggplot2包行可视化处理，以P<0.05为差异有统计学意义。

1.5 DEGs对应PPI网络的构建 蛋白-蛋白互作网络可以基于现有数据对DEGs建立功能性蛋白质关联网络，可通过string数据库(https://cn.string-db.org/)实现，并通过Cytoscape3.7.1软件构建PPI网络，筛选关联度高的hub gene。

1.6 部分Hub gene应用价值的初步验证 利用R语言pROC包，以Cytoscape3.7.1软件中degree值为筛选标准，对部分关联度高的hub gene做ROC曲线，初步验证其作为疾病标志物的价值。

2 研究结果

2.1 DEGs筛选 通过对数据的下载及前期处理，健康对照组相对于病变组做差异表达分析后一共得到18 843个DEGs，根据|LogFC|>1、adj.P<0.01条件筛选得到259个FC值明显改变的基因，包括135个上调基因，124个下调基因(对照组与ICM组)，R语言绘制可视化火山图，见图1。以|LogFC|>1.5、adj.P<0.01为条件共筛选出52个DEGs，包括25个上调基因，27个下调基因，R语言绘制可视化热图，见图2。按|LogFC|的大小分别筛选上调、下调基因前10的DEGs并列表。见表1。

表1 上调、下调DEGs列表

图1 健康心肌组织与缺血性心肌病心肌组织DEGs火山图

图2 健康心肌组织与缺血性心肌病心肌组织DEGs热图(|LogFC|>1.5,adj.P<0.01)

(|LogFC|>1,adj.P<0.01),红色点表示上调基因，蓝色点表示下调基因，黑色点为未达到筛选条件的差异基因。

2.2 DEGs GO分析结果 对筛选得到的135个上调基因，124个下调基因分别做GO富集分析，设定P值<0.05。见图3。上调基因的BP主要富集在炎症反应调节、白细胞游走、创伤修复、细胞之间黏附调节、对内毒素的应答、对细菌来源分子的应答、骨髓白细胞游走、白细胞趋化性、炎症反应的正向调节、单核细胞趋化性；下调基因的BP主要富集在肌系统进程、肌肉收缩、心脏传导系统动态变化。见图3A。上调基因CC主要富集在含胶原纤维的细胞外基质、细胞焦点粘连、细胞基质连接、分泌囊内腔、细胞质囊、囊泡腔、质膜外侧面、初级溶酶体、嗜苯胺蓝粒(淋巴细胞)、血小板α颗粒；下调基因CC主要富集在平滑肌纤维。见图3B。上调基因MF主要富集在糖基化终产物受体(RAGE)受体结合，下调基因MF主要富集在氧化还原酶活性、NADP结合途径、细胞外基质中具有抗压功能的结构成分。见图3C。

2.3 DEGs KEGG分析结果 对筛选得到的135个上调基因，124个下调基因分别做KEGG富集分析，设定P值<0.05。见图4。上调的DEGs主要富集在ECM受体交互作用途径、沙门氏菌感染、HIF-1信号通路、TNF信号通路、Malaria通路、细胞中吞噬体、造血细胞谱系、致病性大肠杆菌感染、细胞的焦点粘连、癌症中的MicroRNAs。下调的DEGs主要富集在药物代谢细胞色素P450途径。

4A 上调基因KEGG富集分析 4B 下调基因KEGG富集分析

2.4 PPI网络的构建 将259个DEGs导入string(https://cn.string-db.org/)数据库构建相互之间的作用网络，在获得的网络中，每个节点(node)代表DEGs表达对应的蛋白质，节点之间的连线(edge)代表蛋白质之间的相互作用，利用NetworkAnalyzer计算无向网络拓扑参数，得出每个节点的连线数量，即degree值，按degree值从大到小筛选节点，degree值越大的节点越有可能是疾病相关的生物标志物，利用CytoHubba插件生成DEGs前10个节点可视化网络图,见图5，并按Degree值大小分别列出前十个节点对应的DEGs、差异上调基因及差异下调基因。见表2。

图5 利用cytoscape构建的DEGsPPI网络图

表2 Degree数值前10位DEGs列表

2.5 部分Hub gene应用价值的验证 选择上述上调、下调基因中degree值最大的前2个hub gene(IL6和RHOA)，结合其在基因芯片中的表达量数值，基于R语言的pROC包绘制ROC曲线图，可见IL6和RHOA的AUC值分别达到了93.3%和96.7%，其作为疾病相关的潜在标志物价值较大。见图6。

6A IL-6 的ROC曲线图 6B RHOA的ROC曲线图

3 讨论

ICM患者与健康人基因表达存在差异，利用生物信息学方法，从ICM患者的基因及蛋白质层面进行挖掘分析，试图找到ICM相关生物标志物，以期对ICM患者的早期诊断、治疗及预后提供参考。与正常对照组相比，DEGs的GO富集分析显示，上调基因的BP富集结果包括肌肉收缩、心脏传导系统的变化，CC富集结果包括平滑肌纤维，MF富集结果包括具有抗压性的细胞外基质成分；下调基因的CC富集结果包括含胶原纤维的细胞外基质(extracellular matrix，ECM)、细胞基质连接、细胞焦点粘连等，下调基因的KEGG富集分析主要富集在ECM受体交互作用途径。ICM患者在支配平滑肌纤维、肌肉收缩、心脏传导等基因中表达上调，这符合心肌缺血病理生理。此外,ICM患者细胞外基质与正常人差异明显，Jiang等[5]研究ICM患者的心外膜脂肪组织(epicardialadpose tissue，EAT)发现，EAT发生了严重的ECM重构，而EAT与心肌是无筋膜隔开的，其与心肌有共同的微循环[6]，可见ECM的变化对ICM的诊治有着重要的意义。此外，ICM患者上调基因中GO富集分析MF富集在氧化还原酶活性、NADP联合，提示ICM患者的发病可能与某些氧化还原酶活性相关，Watanabe等[7]研究发现血浆黄嘌呤氧化还原酶(xanthine oxidoreductase，XOR)活性是冠状动脉痉挛(coronaryspasm ，CAS)发生的独立危险因素，且较高的XOR活性与痉挛的严重程度相关，黄嘌呤氧化酶、黄嘌呤脱氢酶及大部分氧化还原酶均以NAD+/NADH作为电子转移的受体和供体，可见ICM患者的部分氧化还原酶可能成为疾病的潜在生物标志物。

RHOA在ICM患者中表达上调，RHOA表达为蛋白质RhoA，其下游靶点是Rho激酶ROCK(rho-associated protein kinase，ROCK)。在心血管系统中，血管张力的调节主要依赖肌球蛋白轻链激酶活化的Ca依赖性血管收缩和RhoA/ROCK非Ca依赖性血管收缩[8-9]两条途径，RhoA/ROCK可作用于肌球蛋白轻链，使其磷酸化并导致血管收缩，也可直接作用于肌球蛋白轻链磷酸酶，使肌球蛋白轻链去磷酸化受限导致血管收缩[10]。因此RHOA过表达有可能是心肌缺血的重要原因，目前Rho激酶抑制剂盐酸法舒地尔已被广泛应用于临床改善脑组织微循环，包括Y-27643等Rho激酶抑制剂[11]。也有研究[12]表明在进行冠状动脉搭桥手术时，用异氟烷预处理会抑制RhoA/ROCK，激活内皮细胞NO合成酶，导致移植血管冠状动脉前降支血管舒张，对冠状动脉左前降支血管重建后缺血心肌的再灌注非常有用。所以RHOA的过表达可能是ICM心肌缺血的一个潜在靶点，对ICM的预防和早期治疗有意义。

IL-6不仅是一种多效免疫调节细胞因子，可由多种类型细胞分泌[13]。血管内皮细胞中IL-6可以使血小板衍生生长因子(platelet derived growth factor，PDGF)表达升高，PDGF不仅可以刺激血管平滑肌细胞和成纤维细胞的增殖和迁移，同时PDGF也作为一种趋化因子,对单核细胞和中性粒细胞有趋化作用，诱导其产生炎症因子。IL6也可以增加血管内皮细胞的渗透性，促进血管炎与动脉粥样硬化的发生。内皮细胞可以通过L-精氨酸合成、释放NO，调节冠状动脉血管及邻近血管的血管阻力，在心肌缺血梗死后，IL-6可能激活L精氨酸/NO通路[14]。研究[15]显示IL-6对心脏组织的作用因作用的时间和下游信号而不同，在急性缺血损伤期，IL-6信号转导主要起保护心肌组织作用，当长期持续时，则会抑制肌细胞功能造成组织损伤，主要基于JAK-STAT(janus kinase-signal transducers and activators of transcription)通路完成。Podewski等[16]研究发现IL-6在缺血性扩张型心肌病患者的心肌组织中表达下调，且IL-6及其相关因子通过共同作用于gp130、JAKs及STATs为心肌细胞提供重要生存途径，这与筛选结果是一致的。NPPA、NPPB可转录利钠肽，研究显示其在缺血心肌病后期心衰时可升高，目前B型利钠肽、氨基末端B型利钠肽在心力衰竭诊断及预后评估中被广泛应用。TLR2、TLR4主要表达Toll样受体，与炎症相关，部分GO、KEGG富集分析显示ICM与炎症反应、免疫反应相关。本研究显示IL-6 ROC曲线下面积AUC高达93.3%，对ICM的诊断价值较高，应进一步研究其对ICM的诊断和疗效观察价值。

本研究通过生物信息学方法结合高通量数据集，通过健康心肌组织与缺血心肌组织的DEGs分析，并对DEGs进行富集分析、构建PPI网络寻找hub gene及对部分hub gene做ROC曲线验证其应用价值，发现ICM患者RHOA基因的过表达、细胞外基质的变化、氧化还原酶活性、炎症反应的调节等，可能与ICM的发病密切相关，并且本研究发现RHOA和IL-6是ICM潜在的生物标志物，需要进一步验证其用于ICM的诊疗和疗效观察价值。