基于生物信息学探讨合并心力衰竭的扩张型心肌病靶基因的预测

李 果，陈 倩，李鸿炜，关 畅，陈志腾，张玉玲，王景峰

（中山大学孙逸仙纪念医院心内科，广东广州 510120）

扩张型心肌病（dilated cardiomyopathy，DCM）是一种以一个或两个心室增大或扩张，以及收缩力受损为特征的心肌疾病，左心室射血分数（left ventricular ejection fraction，LVEF）通常小于40%，存在收缩功能障碍，可能合并或不合并心力衰竭（heart failure，HF）症状［1］。然而扩张型心肌病是进行性心力衰竭的常见原因，也心脏移植的最常见指征［2-3］，猝死可发生于疾病的各个阶段［4］。扩张型心肌病的发病年龄多为20～60 岁间，在男性中比在女性中更常见［1］，在普通人群中每2 500 人就有1 人发生［5］。近些年来的大量研究表明，扩张型心肌病的最常见病因是特发性的，可能具有家族性或遗传易感性，与细胞骨架，核膜C，收缩蛋白，桥粒，肌节蛋白，线粒体，核膜和RNA 结合蛋白的基因突变均相关，也可能源于各种心肌损伤［1，5］。根据文献报道，扩张型心肌病最常见的遗传易感性就是肌联蛋白（titin，TTN）的截短变异，TTN 编码心肌细胞中表达的巨大的肌联蛋白［6］，可以充当肌节内的生物弹簧，在肌节的形成，机械传感及信号转导中发挥重要作用［7］，其不同剪接形式与心室顺应性有关［8］。此外，参与肌肉收缩功能的编码基因肌球蛋白重链7（myosin heavy chain 7，MYH7），以及编码肌钙蛋白T2（troponin T2，TNNT2）基因等的变异也是扩张型心肌病常见突变类型［9］。病毒感染，心肌炎症，以及细胞内及体液免疫系统的激活导致的自身免疫紊乱等也参与了扩张型心肌病的发生发展［10-11］。目前，利尿剂、ACEI/ARB 或β受体阻滞剂等药物和心脏移植已成为其主要治疗方法，此外有病例报道支持作为ACEI/ARB 的替代药物，血管紧张素受体-脑啡肽酶抑制剂在起始剂量100 mg 时治疗合并心力衰竭的扩张型心肌病患者比ACEI 有更好的疗效和安全性，短期改善心力衰竭效果明显［1，12］。而这些却远远不能满足临床需求［1］。而对于合并终末期心力衰竭的扩张型心肌病，心脏移植是目前唯一可行的治疗方法［13］，新的治疗方法和病理生理机制有待更多探索。随着编程语言的升级和生物数据库的完善，生物信息学分析成为越来越内涵丰富的临床疾病研究手段，在核酸和蛋白质水平展开分析，探究生物大分子的结构及功能信息，为我们提供了疾病进展过程中差异基因发挥作用的更多相关通路注释信息。本研究立足于基因表达汇编（gene expression omnibus，GEO）和人类基因和遗传疾病的在线目录（online mendelian inheritance in man，OMIM）两个公共数据库，运用生物信息学分析方法进一步探究合并心力衰竭的扩张型心肌病的差异表达基因及与疾病相关的信号通路、蛋白互作网络（protein-protein interaction，PPI）分析，以期从分子水平预测可能在合并心力衰竭的扩张型心肌病的疾病进展中发挥重要作用的相关基因。

1 材料与方法

1.1 材 料

从GEO 数据库（https：//www.ncbi.nlm.nih.gov/geo/）下载了以GPL570 Affymetrix Human Genome U133 Plus 2.0 Array 为平台，来源于包含有左室射血分数下降的扩张型心肌病患者和正常心脏功能的心肌组织的2 个基因芯片数据集（GSE29819 和GSE21610），均是由德国波鸿大学Herz-& Diabeteszentrum NRW 实验室测序并上传至GEO 数据库的公共芯片数据集。GSE29819 芯片数据集包括致心律失常性右心室心肌病标本12 例（6 例病例，左右心室心肌组织各取1 例样本），合并心力衰竭的扩张型心肌病标本14 例（7 例病例，左右心室心肌组织各取1 例样本）和正常心脏功能标本12 例（6 例病例，左右心室心肌组织各取1 例样本）；GSE21610 芯片数据集包括未经心室辅助装置支持治疗前和经过心室辅助装置支持治疗后的合并心力衰竭的扩张型心肌病的心肌组织标本各21例（来源于相同的21 例病例），未经心室辅助装置支持治疗前和经过心室辅助装置支持治疗后的合并心力衰竭的缺血性心肌病的心肌组织标本各9 例（来源于相同的9例病例），正常心脏功能标本8例。我们该研究的纳入标准为选择其中只经过标准的抗心力衰竭药物治疗而未经特殊治疗或处理（例如心室辅助装置支持治疗等）的所有扩张型心肌病的心肌组织样本，以及所有正常心脏功能的样本。

1.2 数据的处理与差异基因筛选

采用R 语言对上述两个数据集的标本按照不同芯片分别进行主成分分析（principal component analysis，PCA），观察组别间的分布情况。并用GEO2R 在线工具（https：//www.ncbi.nlm.nih.gov/geo/geo2r/）分析上述各个数据集的差异表达基因，设置筛选条件为log2FC ＞1 或log2FC ＜-1，P＜0.05确定差异表达基因。而其中来自于两张芯片的共同表达差异基因中会包含上下调不一致的基因，直接将所有共同表达差异基因进行与合并心力衰竭的扩张型心肌病相关基因及通路的生信分析，会混入假阳性共同表达基因的影响。为了排除这一混杂因素，且为筛选可以作为临床诊断及预后的预测靶点，相较于正常心脏功能的健康人的表达上调基因显然更有临床应用可行性和更多研究价值，因此我们只选择共同差异基因中表达上调的基因进行分析。对两个数据集获得的差异表达基因分别运用R 语言绘制热图和火山图。并运用Venn Diagramm 对两个数据集筛选得到的差异表达上调基因进行取交集，获得与扩张型心肌病相关的表达一致的差异上调基因。

1.3 共同的差异表达基因的GO 富集及KEGG 信号通路富集分析

提取上述共同的差异表达基因，运用DAVID在线数据库（https：//david.ncifcrf.gov/）进行背景为Homo sapiens 的基因本体（gene ontology，GO）功能富集分析和KEGG（kyoto encyclopedia of genes and genomes，KEGG）信号通路富集分析，并以调整后P＜0.05 为阈值筛查差异基因的主要富集功能和通路。

1.4 共同的差异表达基因的蛋白互助网络分析（PPI）

用在线数据库STRING（https：//string-db.org/）对共同的差异表达基因进行蛋白互作网络的分析，并将所得结果导入Cytoscape 软件，进行可视化和关联性分析，运用蛋白质复合物聚类算法（Molecular Complex Detection，MCODE）插件筛选出关键的蛋白表达分子。

1.5 扩张型心肌病相关的核心基因筛选

在OMIM（https：//omim.org/）中搜索关键词“dilated cardiomyopathy”，寻找文献中已经报道的和扩张型心肌病相关的基因，并和MCODE 分析得到的关键蛋白分子用Venn 图取交集，并定位到阐释我们所取到的交集基因在扩张型心肌病进展中发挥作用的已报道文献。

1.6 用GSEA 差异基因筛选

用基因集富集分析（gene set enrichment analysis，GSEA）方法分别对两个芯片数据集进行核心差异基因筛选，以C2 数据集中的KEGG 基因集为预设基因集进行富集分析。以NOM.P＜0.05 为两个芯片中与扩张型心肌病性状相关的富集通路的筛选标准，将对这些富集通路起关键作用的核心基因选出，分别与OMIM 数据库得到的文献报道的扩张型心肌病相关的候选基因用Venn 图再次取交集，并定位到阐释GSEA 与OMIM 交集的共同基因与扩张型心肌病存在关联的已报道文献。

2 结果

2.1 扩张型心肌病相关的差异表达基因的筛选

经纳入条件的筛选，我们的研究包含来源于芯片数据集GSE29819 的合并心力衰竭的扩张型心肌病标本14 例［病例年龄（57.6±2.3）岁，年龄中位数为57岁，最小年龄为46岁，最大年龄为66岁，3 例男性病例和4 例女性病例，均无冠脉疾病征象，NYHA IV］和正常心脏功能标本12 例［病例年龄（45.2±8.4）岁，年龄中位数为50 岁，最小年龄为23 岁，最大年龄为64 岁，3 例男性病例和2 例女性病例，其中1 例病例年龄和性别均缺失，均因技术原因无法作为心脏移植供体］，来源于芯片数据集GSE21610 的未经心室辅助装置支持治疗的合并心力衰竭的扩张型心肌病标本21 例［病例年龄（48.6±2.9）岁，年龄中位数为50 岁，最小年龄为16 岁，最大年龄为66 岁，19 例男性病例和2 例女性病例，均无冠脉梗塞征象，LVEF 均小于40%］和正常心脏功能标本8 例［病例年龄（29.0±6.2）岁，年龄中位数为25.5 岁，最小年龄为5 岁，最大年龄为64 岁，6 例男性病例和2 例女性病例，均无心脏疾病历史且均因技术原因无法作为心脏移植供体］。

这两张芯片分别进行的PCA 主成分分析用散点图展示，散点图中每一个点分别代表一个样本，而正常心脏组合扩张型心肌病组在两个数据集中均显示较为一致的聚类结果（图1A、B）。根据筛选条件，在火山图中将筛选得到的有意义的差异表达基因用不同颜色进行展示（图1C、D）。

图1 芯片数据集PCA 图和差异表达基因火山图Fig.1 PCA plots of gene chips and Volcano maps of differentiate expression genes

这些结果包含我们重点关注的在扩张型心肌病组表达上调的差异基因，GSE29819 得到表达上调的差异基因843 个，GSE21610 得到表达上调的差异基因774 个。这些差异表达上调的基因在各个样本中的表达通过热图进行展示，可以看到各组样本的差异基因表达较为一致（图2A、B）。通过Venn Diagram 将两个数据集的差异表达上调基因取交集，得到173 个与扩张型心肌病相关的共同表达上调的差异基因（图2C）。

图2 表达上调差异基因热图和韦恩图Fig.2 Heatmaps and Venn diagram of up-regulated genes

2.2 扩张型心肌病相关的共同差异表达基因的GO 及KEGG 信号通路富集分析

运用DAVID 在线数据库对筛选得到的共同差异表达基因进行背景为Homo sapiens 的富集分析，得到关于GO 的富集信息。表达上调的差异基因富集到的“生物过程”主要包括炎症信号和细胞增殖分化的正负调控（图3A，表1）；在“细胞组成”方面，差异基因主要富集在细胞外基质等区域（图3B，表2）；而“分子功能”方面，则提示差异基因主要与Wnt 蛋白、转化因子受体等生物大分子的结合及生长因子及Wnt 受体活化相关（图3C，表3）。在KEGG 信号通路的富集分析中显示，表达上调的差异基因主要与Wnt，FoxO 等经典的凋亡信号传导通路相关（图3D，表4）。

2.3 共同差异表达基因的蛋白互助网络构建及模块分析

用String 对共同差异表达基因进行蛋白互作网络构建。并将得到的结果导入Cytoscape 软件对蛋白互作网络进行可视化分析与筛选。利用Cytoscape 中的Network Analyser 工具对蛋白互作网络中的各个节点进行无方向性的分数计算，得到各个节点的Degree 值。以节点大小表示Degree的值，节点颜色从红到绿表示各节点Neighborhood Connectivity 的从高到低，边的粗细表示边的combined_score 值，并用Attribute Circle Layout 的布局对所有的蛋白节点进行排布，并将Degree≥4的节点围在了内层（图4A）。设置筛选条件为Degree≥4 对蛋白节点进行过滤，得到31 个重要的蛋白分子。并利用MCODE 插件对这31 个重要的蛋白分子，以默认参数Node Score Cutoff 为0.2，K-core 为2，Max.Depth 为100 进行聚类关联分析，得到3 个得分较高的Cluster，分别是Cluster1（6 个节点，得分6，图4C），Cluster 2（4 个节点，得分4，图4D），Cluster 3（9 个节点，得分3.75，图4E）。

表1 GO 生物过程富集分析结果Table 1 Results of GO Biological Process analysis

表2 GO 细胞组成富集分析结果Table 2 Results of GO Cell Component analysis

表3 GO 分子功能富集分析结果Table 3 Results of GO Molecular Function analysis

表4 KEGG 分子功能富集分析结果Table 4 Results of KEGG analysis

图3 表达上调差异基因GO 和KEGG 富集分析结果Fig.3 GO and KEGG analysis results of up-regulated genes

2.4 扩张型心肌病相关的核心基因筛选

在OMIM 人类基因和遗传疾病的在线目录中搜索到文献中已报道的与扩张型心肌病相关的基因757 个，并和MCODE 分析得到的19 个关键蛋白分子用Venn 图取交集，得到3 个文献中已报道的候选基因，分别为NPPA（natriuretic peptide precursor A），ApoA1（apolipoprotein A1）和COL6A1（collagen VI，alpha-1 polypeptide）（图4B）。

2.5 用GSEA 差异基因筛选

图4 蛋白互作网络和候选基因韦恩图Fig.4 Protein-protein interaction network and Venn diagram of candidate genes

GSEA基因集富集分析以C2数据集中的KEGG基因集为预设基因集，分别对两个芯片数据集进行核心差异基因筛选，并以NOM.P＜0.05 为与扩张型心肌病性状相关的富集通路的筛选标准，分别在GSE21610 芯片中得到10 条KEGG 通路，GSE29819 芯片中得到3 条KEGG 通路，并将这些通路用散点图进行展示，横纵坐标分别是通路的标准化P 值和KEGG 通路名称，不同性状分别代表了来源于不同芯片数据集，点的大小代表了富集到的基因数量，点的颜色代表了ES 值，及基因在该通路中的富集情况，ES 值越高说明这些基因在通路中富集越好，而不是散在分布，可以说明这些基因在该通路中有共同的表达趋势（图5A）。并且可以看到两个芯片共同富集到了β-丙氨酸代谢通路。

将这些富集通路中的核心基因分别于OMIM数据库得到的文献中报道的扩张型心肌病相关基因取交集，Venn 图中两个基因集与OMIM 共同存在的基因2 个，分别为PRKCA（protein kinase C，alpha polypeptide）和BMP2（bone morphogenetic protein-2，图5B）。

在GSE29819 芯片中，PRKCA是紧密连接通路的核心基因，BMP2是Hedgehog 信号通路的核心基因（图5C）；在GSE21610 芯片中，PRKCA是Wnt信号通路，黑色素生成和非小细胞肺癌的核心基因，BMP2是基底细胞癌的核心基因（图5D）。

3 讨论

扩张型心肌病是以收缩功能障碍为特征的具有遗传易感性的心肌疾病，伴或不伴心力衰竭症状，其发生发展机制与自身免疫紊乱及肌节或细胞骨架蛋白基因的突变等相关［1，11］。其中编码巨大肌联蛋白的TTN 基因的截短变异［6］，以及参与肌肉收缩功能的肌球蛋白重链MYH7，和肌钙蛋白TNNT2等基因的变异均是扩张型心肌病的常见遗传易感类型［9］，与扩张型心肌病的致病机理密切相关。然而合并心力衰竭的扩张型心肌病是临床治疗的难题，除了心脏移植外尚无可靠的治疗方案来逆转心脏功能。

图5 GSEA 的KEGG 预设基因集富集结果Fig.5 GSEA results by KEGG gene sets

在本研究中，我们立足于对GEO 数据库的分析，选取了以GPL570 为研究平台的两张基因数据集芯片，以合并心力衰竭的扩张型心肌病的心肌活检样本为实验组，正常心脏功能的心肌活检样本为对照组，将两张芯片中运用GEO2R 在线差异基因分析方法筛选与合并心力衰竭的扩张型心肌病相关的表达上调的差异基因通过Venn 图取交集，得到173 个表达上调的差异基因为候选基因。对候选基因的GO 和KEGG 分析主要集中在炎症信号和细胞增殖分化的正负调控，细胞外基质及Wnt 蛋白、转化因子受体等生物大分子的结合及生长因子及Wnt 受体活化等相关生物学过程，并且主要与Wnt，FoxO 等经典的凋亡信号传导通路相关。

同时将173 个候选基因导入STRING 进行蛋白互作网络分析，并进行MCODE 模块评分分析，得到了19 个可能与合并心力衰竭的扩张型心肌病的发生发展相关的基因。而这些基因属于胶原蛋白（collagen，COL）家族（COL8A1，COL14A1，COL6A1），钠尿肽（natriuretic peptide，NPP）家族（NPPA，NPPB，NPR3），F-BOX 家 族（FBXO40，FBXO32，FBXL22），与骨相关的基因：骨膜蛋白（periostin，POSTN），骨形态发生蛋白4（bone morphogenetic protein 4，BMP4），BMP4 拮抗剂（chordin like 1，CHRDL1），以及载脂蛋白A1（apolipoprotein A1，ApoA1），釉质素（enamelin，ENAM），脑啡肽原（proenkephalin，PENK），胰岛素样生长因子结合蛋白3（insulin-like growth factor-binding protein 3，IGFBP3），锌指蛋白213（ring finger protein 213，RNF213），皮屑蛋白样蛋白1（leprecan like protein 1，LEPREL1）和膜金属肽链内切酶（membrane metalloendopeptidase，MME）。

根据文献报道，Al-Yacoub 等［14］的实验证明F-Box 蛋白家族的成员FBOXO32是引发扩张型心肌病的新型突变基因，FBOXO32的错义突变影响高度保守的氨基酸，进而被认为会严重损害与SCF 蛋白的结合，同时该实验还证明了有FBXO32突变的患者的心脏存在调节自噬的蛋白的积累，提示着FBOX32是扩张型心肌病发生发展中的重要一环，可能成为扩张型心肌病治疗的潜在靶基因。由于Ⅵ型胶原蛋白基因的突变会导致肌肉发育不良，Crossman 等［15］的实验通过Western-Blots鉴定特发性扩张型心肌病的心肌组织中VI 型胶原蛋白表达与正常心肌相比增加了2.4 倍，和横向小管的肥大。研究表明，横向小管是质膜的延伸部分，可以与肌浆网在心肌细胞表面形成与钙释放耦联心肌收缩相关的电连接，其重塑及肌营养不良蛋白变化导致的纤维化均与晚期心力衰竭相关［16-17］。而ⅩⅢ和XIV 类型的胶原纤维与扩张型心肌病的关系尚未有文献报道，根据我们生物信息学预测的结果表明COL8A1，COL14A1在合并心力衰竭的扩张型心肌病中均有明显升高，很可能提示着他们和COL6A1一样在扩张型心肌病的心力衰竭发生发展中发挥重要作用，有望成为扩张型心肌病失代偿期的预后评估指标。根据文献报道，合并心力衰竭的扩张型心肌病患者会产生低水平的心脏选择性丝氨酸蛋白酶corin［18］，而corin参与了钠尿肽前体的裂解和活性形式的生成［19］，因此corin 和钠尿肽家族的成员均是临床上心力衰竭的常用的心脏损伤标志物，可以用来提示扩张型心肌病失代偿期的发生，及心力衰竭症状的严重程度［18］，与我们的样本选择和分组及生物信息学预测结果相一致。

由于传统筛选表达差异基因的方法会掩盖一些具有重要生物学意义但表达上调倍数较低而被掩盖的基因，我们又采用了GSEA 富集分析的方法来筛选与合并心力衰竭的扩张型心肌病相关的KEGG 信号通路及其核心基因，以NOM.P＜0.05为筛选标准，从两张芯片中我们分别得到10 条和3 条KEGG 信号通路，主要集中在Wnt，JAK-STAT等信号通路及β-丙氨酸代谢，核心基因158 个和46 个。此外，我们还通过OMIM 人类基因和遗传疾病数据库筛选出了文献中已报道的与扩张型心肌病相关的基因。最终，我们将GEO2R 和GSEA得到的候选基因分别与OMIM 库中的已报道基因用Venn 图取交集，分别得到了NPPA，ApoA1和COL6A1三个表达量上调的基因，PRKCA和BMP2两个与合并心力衰竭的扩张型心肌病相关的KEGG 信号通路的核心基因。Hu 等［20］的研究表明，编码蛋白激酶Cα（PRKCA）的基因有增强心肌收缩力的功能，可能与扩张型心肌病的心脏功能失代偿期心肌收缩力的增加有关。Ye 等［21］证明，骨形态发生蛋白2（BMP2）可能通过Kielin/chordin样蛋白来促进扩张型心肌病晚期心力衰竭的发生发展。

然而Sezgin 等［22］的实验表明，ApoA1在不稳定性扩张型心肌病和稳定性扩张型心肌病中与正常心肌组织中相比显著降低，在不稳定性扩张型心肌病中降低尤为明显。文献提示ApoA1的降低与扩张型心肌病病人血清中的细胞因子IL-6 等的升高有关，提示疾病的不良预后，而这与我们的生信分析得到的结果相反，这可能与我们纳入的合并心力衰竭的扩张型心肌病病例都经过延缓扩张型心肌病病程进展的标准药物治疗，体内与炎症相关的细胞因子表达降低，因此ApoA1表达相应地升高相关。这进一步提示了ApoA1在扩张型心肌病的病程进展的不同时期可能存在不一致的表达水平，其在疾病进程中发挥的作用有待更深入的研究。

由于生物信息学对基因芯片分析方法的局限性，我们只能对ApoA1在扩张型心肌病合并心力衰竭的疾病发生发展过程中出现的表达水平不一致的原因提出合理猜想，ApoA1与细胞因子相互作用的具体机制及其在疾病进展中参与的信号通路需要更多的探究与实验。同样，我们的实验发现的PRKCA，BMP2，NPPA，COL6A1等基因与扩张型心肌病合并心力衰竭的疾病发展相关，并通过KEGG 及GO 等方法预测了其可能参与的信号通路，但是其在临床疾病发生过程中的具体作用靶点及相关线索仍需进一步的研究。

综上所述，F-Box 家族，胶原蛋白家族的基因，以及蛋白激酶Cα和骨形态发生蛋白2 等都在扩张型心肌病合并心力衰竭的疾病进程中从心肌自噬蛋白积累，肌横向小管重塑及肌营养不良蛋白增加，和增强心肌收缩力等方面发挥重要作用，而钠尿肽家族的基因可以用来提示扩张型心肌病失代偿期的病情严重程度及心肌损伤程度。根据生物信息学分析筛选出的候选基因PRKCA，BMP2，NPPA，COL6A1等，可通过心肌收缩相关调节或自噬蛋白积累等方式对扩张型心肌病合并心力衰竭期的疾病进程产生影响，其可能在晚期心功能失代偿的临床疾病进展中发挥着重要作用，为临床预后的判断及治疗提供有意义的研究线索及方向。