冠心病差异表达基因的鉴定及生物信息学分析

李召水,王光静,乔友进,生伟,黄强,池一凡

(1 青岛大学附属青岛市海慈医院心外科,山东 青岛 266033；2 青岛大学附属青岛市市立医院心外科)

心血管疾病是目前世界上人类死亡的主要原因之一。冠心病(CHD)是最常见的一种心血管疾病，在全球范围内每年导致超过700万人死亡[1]。2014年的一项研究显示，近1/5的男性和1/10的女性死于CHD[2-4]。据估计，未来20年CHD患病率将增加约10%[5]。CHD已成为威胁人类健康的重要疾病之一，对CHD发病机制及有效疗法的研究和探索从未停止。CHD的主要危险因素包括血脂异常、糖尿病、动脉硬化、肥胖、吸烟、久坐的生活方式、压力、年龄、男性和家族病史等[6]，但其具体发病机制尚不完全清楚。既往研究显示，遗传因素在心血管疾病的发生过程中发挥了极大的作用[7]。CHD发生的生物学机制有多种，其中研究较为清楚的为炎症反应，炎症反应失调是CHD发生的一种潜在的生物学机制[8]。相关研究表明，基因表达差异，尤其是炎症调控相关基因表达异常与CHD的发生紧密关联[9]。除了炎症异常调控外，还有其他因素的变化参与CHD的发生。本研究旨在通过对GEO数据库中CHD发生的基因表达谱进行生物信息学分析，揭示与CHD疾病发生相关的生物学过程及信号通路，为进一步阐明CHD的发病机制提供有价值的信息，并为CHD的诊断、治疗提供新的思路。

1 资料和方法

1.1 基因微阵列数据收集

CHD样本的基因表达芯片来自GEO数据库(http://www.ncbi.nlm.nih.gov/geo/)[10-12]。以“Coronary Heart Disease”为关键词在GEO数据库中进行检索，最终在210个相关数据集中选取2个来自Affymetrix Human Genome U133 Plus 2.0 Array分析平台的基因集GSE71226及GSE19339，共包含7个CHD样本和7个正常样本的表达矩阵。

1.2 原始数据预处理及差异表达基因(DEGs)鉴别

下载原始数据，采用R语言(affy，limma包)对其进行噪声去除、分位数归一化等处理，然后筛选CHD组和正常对照组的DEGs。筛选DEGs的阈值设定为P值<0.05，且|log2(Fold change)|≥1。最后，采用R语言(pheatmap包)对基于mRNA表达水平的组样本进行可视化层次聚类分析。

1.3 GSE71226和GSE19339中共有差异表达基因(co-DEGs)的Venn分析

通过Draw Venn Diagram线上数据库(http://bioinformatics.psb.ugent.be/webtools/Venn/)对GSE71226和GSE19339数据集中的co-DEGs进行分析[13-15]。将待分析基因列表上传到数据库，即可显示维恩图及相关共有基因列表。

1.4 基因本体(GO)和基因通路富集分析

GO注释分析通常用于大规模转录组数据的功能研究。KEGG(Kyoto Encyclopedia of Genes and Genomes)包含了多种生物化学通路。将待分析基因列表上传至DAVID生物信息学资源6.8数据库(https://david.ncifcrf.gov/)[16-17]，即可显示GO及 KEGG分析结果，将其下载为文本文件。最后，通过R语言(ggplot2包)可视化GO结果。

1.5 基因集富集分析(GSEA)

将特定规格的矩阵表格加载到GSEA_4.0.2软件，通过GSEA online进行可视化即可完成GSEA分析[18-19]。DEGs途径富集的阈值为P值<0.01。

1.6 蛋白质调控网络分析

DEGs的蛋白质相互作用(PPI)网络分析通过STRING (https://string-db.org/)在线分析软件完成[19-20]。将基因列表上传到多个蛋白质分析菜单栏，稍后即可显示PPI结果。最后用Cytoscape软件将具体的网络图可视化。

1.7 实时荧光定量聚合酶链反应(RT-qPCR)验证CHD病人核心DEGs的表达水平

取青岛市市立医院心脏外科10例50～80岁CHD病人和10例同龄健康人的外周血样本，使用高效血液总RNA提取试剂盒(天根生化科技(北京)有限公司，Lot#DP443)提取总RNA。用Oligo(dT)引物(Takara,cat#3806，Lot#T2301AA)在65 ℃条件下退火5 min得到mRNA，用RevertAid逆转录酶(Thermo Scientific，#EP0441)和dNTP混合物(Takara，Cat#4019,Lot#AI11312A)进行逆转录得到cDNA模板。最后使用PowerTrack SYBR Green Master Mix(Thermo Scientific，#4367659)及基因特异性引物进行RT-qPCR，检测目的基因的相对mRNA水平。引物序列见表1。

表1 RT-qPCR引物序列

2 结 果

2.1 CHD DEGs的筛选

从GEO数据库中收集了7例CHD病人和7例正常对照者的mRNA表达谱。根据|log2(Foldchange)|≥1、P值<0.05的筛选条件，GSE71226数据集中共鉴定出2 262个DEGs，其中包含上调基因694个及下调基因1 568个(图1A)；GSE19339数据集中共鉴定出537个DEGs，其中包含上调基因263个及下调基因274个(图1B)。对这些DEGs进行热图聚类分析结果显示，CHD组和正常对照组的基因表达模式差异显著(图1C、D)。

由于样本来源不同(GSE71226数据集中样本来自CHD病人和正常人的外周血；GSE19339数据集中样本分别来自经皮冠状动脉介入治疗的CHD病人冠状动脉闭塞部位的血管和正常人外周血)，两个数据集中分析得到的DEGs具有一定差别。而且，两个数据集中病人信息极少，故无法分析年龄、性别和病史对CHD DEGs的影响。

A、B分别为火山图显示GSE71226和GSE19339数据集中CHD的DEGs。红色的点代表表达上调基因，绿色的点代表表达下调基因；垂直黑线分别为上下2.0倍，水平黑线表示P值为0.05(双尾不配对t检验)。C、D分别为热图显示GSE71226和GSE19339数据集中CHD病人与正常人的DEGs分布。红色条码代表CHD病人表达上调的基因，蓝色条码代表CHD病人表达下调的基因；右侧统计图为上调基因(上)和下调基因(下)的平均相对表达水平，经t检验显示差异有统计学意义。

2.2 co-DEGs的Venn分析

为了较为精确地研究CHD的DEGs，本研究分析了GSE71226和GSE19339两个数据集中的co-DEGs。结果筛选得到两个数据集中共同上调基因8个及共同下调基因114个，共计122个co-DEGs。见图2。两个数据集中大部分co-DEGs均为表达下调基因，提示这些共同下调基因可能是CHD发病的关键基因。

GSE71226和GSE19339数据集中共同上调基因(左)与共同下调基因(右)的Venn重叠图，二者交集中的数字即为co-DEGs数目。UG：上调基因；DG：下调基因。

2.3 CHD co-DEGs的GO和KEGG分析

为了阐明DEGs的生物学功能，对以上122个co-DEGs进行了GO富集分析。结果显示，CHD中大多数的co-DEGs参与的生物学过程(biological process)为mRNA加工和剪接调控、细胞内转录调控(图3A)；co-DEGs所属的细胞成分(cell components)为核质、细胞核和细胞质(图3B)；其分子功能(molecular functions)主要为poly(A)RNA结合、蛋白结合、DNA结合(图3C)。KEGG富集分析显示，大多数co-DEGs显著富集的信号通路为剪接体(图3D)。以上分析结果表明，CHD的发生与细胞整体蛋白质表达调控紊乱或RNA剪接紊乱具有重要关联。

2.4 CHD DEGs的GSEA分析

为了进一步分析CHD DEGs可能参与的信号通路，本研究对其进行了GSEA分析。结果显示，两个GEO数据集中DEGs共同低表达的基因富集的信号通路为mRNA过程的调节及DNA损伤修复(DNA damage repair)(图4)。表明CHD发生过程中，涉及mRNA调节过程及DNA损伤修复途径的相关基因表达水平下降。GSEA分析结果与GO分析结果相一致。

2.5 CHD DEGs的蛋白质调控网络分析

为了筛选CHD的关键DEGs，本研究对122个co-DEGs进行了PPI分析。结果显示，两个数据集共同下调的基因大部分处于PPI网络中间，而共同上调的基因则处于网络边缘。其中，位于PPI网络中心的基因分别为LUC7L3、HNRNPA1、SF3B1、ARGLU1、SRSF5、SRSF11、SREK1、PNISR、DIDO1、ZRSR2及NKTR(图5)。提示这些基因的异常低表达可能在CHD的发生过程中发挥了重要作用。

A、B、C分别为GO分析的biological process、cell components和molecular function；D为KEGG分析。条形图右边的数字表示每一项的基因数量。

GSE19339(A)和GSE71226(B)数据集中共同低表达的基因富集的信号通路。

2.6 CHD关键DEGs分析

筛选出的11个关键DEGs在CHD中均显著低表达(图6)。GO分析结果显示，这些关键DEGs涉及的生物学过程为DNA转录和RNA剪接体调控(表2)。表明CHD的发生与RNA剪接异常调控具有重要关系。

表2 CHD关键DEGs参与的生物学过程

2.7 RT-qPCR验证CHD关键DEGs的表达

分别收集10例CHD病人及10例正常人的外周血，对筛选出的关键DEGs的表达水平进行了RT-qPCR验证。结果显示，CHD病人外周血中这些DEGs的表达水平均较正常人显著下调。见表3。

表3 RT-qPCR验证CHD关键DEGs的表达

3 讨 论

尽管对CHD进行了40多年的基础和临床研究，但其具体发病机制仍不完全清楚。通过分析CHD发生过程中涉及的生物学途径，增加对CHD发病机制的了解，可为CHD的临床治疗及预后判断提供新思路。

剪接体被证明是一种蛋白质定向金属酶[21]。作为真核细胞中最复杂的调控机制之一，剪接体从初级转录本中去除内含子序列，生成功能性mRNA和长链非编码RNA(lncRNA)[22]，这一过程称为选择性剪接。选择性剪接是一个动态且受调控的生物学过程，受到一系列变量的影响，如顺式调控序列和反式作用因子、转录过程和DNA/RNA的甲基化等[23-24]。多项研究表明，异常可变剪接与人类疾病有关，它既可能是疾病的发生原因，也可能是疾病造成的结果[25]。有研究结果表明，参与剪接体正常功能的基因突变被认为是脊髓性肌萎缩、色素性视网膜炎和普瑞德-威利综合征等的关键因素[26-28]。然而，剪接因子中导致人类心脏病变的突变并不多见。到目前为止，只有剪接因子RNA结合基序蛋白20(RBM20)的突变被证实与心脏病有因果关系[29-31]。此外，相关研究结果表明，与RNA剪接相关基因在心脏病中异常表达。例如，剪切因子SF3B1在患病的人和小鼠心脏中均表达上调[32]，Rbfox1基因在人类和小鼠心脏中表达下调[33]。然而，CHD病人中DEGs一直未被明确阐述。

本研究分析了GEO数据库的GSE71226和GSE19339数据集中CHD病人的基因表达数据，拟筛选与CHD发生密切相关的DEGs，探讨CHD基因水平的发病机制。结果显示，1.118%～2.954%(GSE71226：2.954%；GSE19339：1.118%)的基因表达水平上调，同时有1.165%～6.667%(GSE71226：6.667%；GSE19339：1.165%)的基因表达水平下调，表明CHD的发生与细胞中基因表达的变化密切相关。由于样本来源和各微阵列平台研究都存在差别，综合分析各种微阵列数据集可以获得更为准确的结果，故选择了两个数据集中8个共同表达上调基因及114个共同表达下调基因进行进一步分析。GO注释分析结果表明，这些DEGs参与了DNA转录和mRNA剪接调控，提示CHD的发生与细胞中RNA剪接紊乱有关。选择性剪接是一种可实质上改变基因表达模式的转录后机制。高达95%的人类基因具有多外显子可变剪接形式，表明可变剪接是人类基因组功能复杂性的最重要组成部分之一。本研究结果表明，大部分的CHD DEGs是可变剪接相关的基因，提示可变剪接调控在心脏病的研究中应受到更多的重视。

红圈表示上调的基因，蓝圈表示下调的基因。

GSE71226(A)和GSE19339(B)数据集中关键DEGs相对于正常人的表达水平，P值显示在相应的条形图上面。

在DEGs调控网络中，表达下调的LUC7L3、HNRNPA1、SF3B1、ARGLU1、SRSF5、SRSF11、SREK1、PNISR、DIDO1、ZRSR2和NKTR位于网络控制中心，且均为DNA转录和RNA剪接调控相关基因。既往研究表明，LUC7L3通过RE和RS域参与了剪接体的形成，在心脏钠通道剪接调节人类心力衰竭中发挥作用[34-35]。HNRNPA1为异质性核糖核蛋白(hnRNP)复合体中含量最丰富的核心蛋白之一，在选择性剪接的调控中发挥关键作用。SF3B1为一种重要的pre-mRNA剪接因子，与癌症突变相关，并可以作为靶向药物治疗靶点[36-40]。在剪接体装配的早期阶段，SF3B1在pre-mRNA剪接位点的小核核糖核酸蛋白(snRNP)之间促发了一系列依赖ATP的结构和成分重排，最终完成pre-mRNA剪接的行为[36,41-42]，但其在CHD中的作用尚未得到证实。据报道，ARGLU1为一种转录共激活因子和剪接调节因子，对应激性激素信号转导和发育以及多种癌症调控非常重要[43-44]。SRSF5是pre-mRNA剪接因子中SR的家族成员，是剪接体的一部分[45]。已有研究结果表明，SRSF5作为一种新型的致癌剪接因子，在多种癌症和免疫调节中发挥关键作用[46-51]，但其在CHD中的作用未见报道。SRSF11为一种在可变剪接过程中发挥作用的剪接因子[52]。SREK1为富含SR剪接蛋白家族的一个成员[53]。PNISR，又被称为SFRS18，使用公开交互数据库的数据挖掘也支持了LUC7L3和SFRS18在RNA剪接中的相互作用[54]。GARCIA-DOMINGO等[55-56]研究表明，DIDO1通过上调procaspase 3和9参与细胞凋亡的激活。此外，FÜTTERER等[57]观察到，小鼠中DIDO的缺失与骨髓增生异常综合征相关。FLEISCHMAN等[58]的研究则表明，ZRSR2突变病人的常见临床特征为白细胞减少、血小板减少或骨髓母细胞百分比增加的大细胞性贫血。本研究中筛选到的CHD DEGs大部分都是mRNA剪接相关基因，这些基因通过RNA剪接功能调控不同的人类疾病。但是，这些基因与CHD之间的关系目前尚未被报道。

综上所述，本文结果显示，CHD病人RNA剪接相关基因的表达水平发生显著改变，表明RNA剪接调控在CHD的发生过程中可能发挥了重要作用，但其在CHD中的具体作用机制仍有待进一步研究。本研究结果为CHD的进一步研究及高危人群的筛查提供了新的思路。