黄 蕾，秦显雨，吴岳恒，林吉进

［1.南方医科大学第二临床医学院，广州 510515；2.广东省心血管病研究所心内科广东省人民医院（广东省医学科学院），广州 510080］

肥厚型心肌病（hypertrophic cardiomyopathy，HCM）是最常见的由常染色体显性遗传的心肌病，既往普遍认为其发病率在普通人群中为1∶500～1∶200［1-2］。HCM 是年轻人包括训练有素的运动员猝死的最常见原因之一［3］。因此，HCM 发病机制的研究对HCM 的诊疗有着极其重要的意义。有研究表明，能量代谢、血管形成和炎症反应在心肌病理性肥大和心肌重构过程中起着至关重要的作用［4-5］。但其具体机制至今尚未明确，系统性分析其分子机制及转录图谱有利于进一步揭示疾病的整体病理进程。本研究拟通过系统性分析HCM患者的高通量测序芯片，构建心肌病的差异性基因表达谱，进行生物信息学注释，从而发现与HCM 相关的关键基因及信号通路，并进一步阐明其潜在的分子机制。

1 材料和方法

1.1 材 料

本研究所使用的高通量测序数据集（GSE36961）来自于美国国家生物信息中心（NCBI）下的高通量基因表达数据库Gene Expression Omnibus（https：//www.ncbi.nlm.nih.gov/geo/），平台编号为GPL15389（Illumina HumanHT-12 V3.0 expression beadchip），表达数据类型为Expression profiling by array，种属为Homo sapiens。该数据集中包含106例HCM 样本及39 例对照组样本，HCM 的基因组样本来自于HCM 患者的外科手术切除组织，对照组样本来自于正常人的心肌组织。

1.2 方 法

1.2.1 差异表达基因谱构建 通过使用R 语言中Bioconductor 工具集中的“limma 程序包”，来完成GSE36961 数据集中的HCM 和正常心肌组织间的差异表达基因分析，并应用贝叶斯算法对P值进行校正，定义|LogFC|≠0 且校正后P＜0.05 为标准来对数据集进行筛选，从而构建HCM 患者的差异表达基因谱。

1.2.2 加权重共表达网络构建及模块识别 通过使用R 语言中的“WGCNA 程序包”，来进行加权重共表达网络的构建及模块的识别。首先剔除离群样本从而使网络构建结果稳定，然后选择适当的软阈值β使构建的网络近似为一个无尺度网络，并将数据转化为拓扑重叠矩阵。随后利用动态剪切树法行层次聚类，从而得到差异基因的整体聚类树，聚类树的分支即定义为基因模块。然后计算这些基因模块的特征向量基因（module eigengene，ME），通过聚类合并这些模块的ME，可得到相似程度高的基因模块，从而构建出差异表达基因的加权重共表达网络。接着对模块的ME 进行分析，最终得到我们所需要的目的模块。

1.2.3 模块内基因功能富集分析 为了解识别出来的模块内基因所具有的生物学功能，本研究使用ToppGene 数据库（https：//toppgene.cchmc.org/）进行基因本体论分析（gene ontology，GO）及京都基因与基因组百科全书分析（Kyoto encyclopedia of genes and genomes，KEGG）。通过GO 分析，我们可得到各模块中基因所富集的生物过程（biological process，BP）、细胞组分（cellular components，CC）及分子功能（molecular function，MF）。通过KEGG分析，我们可得到各模块基因所富集的信号通路。以上均定义P＜0.05 为显著富集。

1.2.4 蛋白质相互作用网络构建 使用STRING数据库（https：//string-db.org/）对差异表达基因进行蛋白质相互作用网络分析。在蛋白互作网络中，一个节点（Node）代表一个蛋白质，节点间的连线（Edge）代表蛋白间的相互作用。使用Cytoscape软件（3.7.1 版）对蛋白质互相作用的结果进行绘图，并使用软件中的CytoHubba 插件通过12 种方法筛选出连接度最高的前10 个关键基因。

1.3 统计学分析

本研究中的差异表达基因谱的构建、加权重共表达网络的构建及模块识别均使用基于R 语言平台的R 软件（3.6.1 版）来完成。以上均定义P＜0.05 为差异有统计学意义。

2 结果

2.1 肥厚型心肌病差异表达基因谱构建

通过运用R 语言中的“limma 程序包”以|Log-FC|≠0、P＜0.05 为标准对数据集GSE36961 进行分析，筛选出8 002 个差异表达基因。根据差异表达基因的结果构建火山图（图1），从而对样本之间和基因之间的差异表达结果进行直观展示。同时对前40 个显著性最高的差异表达基因进行层次聚类分析并构建热图（图2），热图的横坐标代表组织样本的聚类，一列代表一个样本，其聚类基于样本间基因表达的相似性，因此样本间基因的表达越接近，则靠得越近；纵坐标则代表差异表达基因的聚类，一行代表一个基因，聚类基于基因在样本中表达的相似性，因此基因在样本中的表达越接近，则靠得越近。由此可见，表达模式相近的基因具有功能相关性。

2.2 加权重共表达网络构建及模块识别

如方法学所示，使用R 语言中的“WGCNA 包”来实现差异表达基因的加权重共表达网络的构建。在选择合适的软阈值β以构建无尺度网络时，需考虑对各基因节点（Node）的平均连接度、中位连接度及最大连接度进行适度保留，本研究最终选取β=11 用于构建加权重共表达网络。确定β值后，通过动态剪切树法进行模块的初步识别并合并ME 相似的模块，去除灰色模块（此模块由无法分配至任何一个模块的基因构成）后，最终得到11 个基因共表达模块。

在纳入研究的临床信息中，性别（gender）和年龄（age）为非时间相关的变量。以ME 及其相应变量的Pearson 相关系数，即模块隶属度（module membership，MM），来代表模块与相应临床信息的相关性，可见青色模块MEcyan（Pearson cor=0.77，P=4e-29）、紫红色模块MEmagenta（Pearson cor=0.76，P=2e-28）均与HCM 中等程度相关，其余模块与临床信息的关系均为弱相关或不相关（图3）。各模块与性别、年龄的相关系数绝对值均小，且相关性均基本为不显著。

为筛选出基因显著水平（gene significance，GS）与MM 有较高相关程度且在模块中发挥更重要作用的基因，将MM 相对于GS 做散点图，可见青色模块中cor＝0.68，P＜1e-200（图4），紫红色模块中cor＝0.76，P=9.8e-73（图5）。

2.3 模块内基因功能富集分析

将两个模块中筛选出来的基因分别导入ToppGene 数据库（https：//toppgene.cchmc.org/）中进行GO 分析及KEGG 通路分析，设置P＜0.05 为显著富集的标准，富集相关性最强的前10 个条目绘图。通过GO 分析，可见青色模块基因的生物过程主要富集在能量代谢中（图6A）；细胞组分主要富集在线粒体中（图6B）；分子功能则主要富集在能量代谢相关蛋白活性及结合中（图6C）。而紫红色模块通过GO 分析，可见其生物过程主要富集在血管形成中（图7A）；细胞组分主要富集在含细胞外基质相关成分中（图7B）；分子功能则主要富集在细胞外基质相关蛋白活性及结合中（图7C）。通过KEGG 分析，可见青色模块中基因的信号通路主要富集于能量代谢中（图6D），紫红色模块基因的信号通路主要富集在细胞外成分形成的相关通路中（图7D）。

2.4 蛋白质相互作用网络的构建

将两个模块中所有筛选出来的基因（青色模块2 816 个，紫红色模块380 个，共3 196 个基因）导入STRING 数据库（https：//string-db.org/），得到差异表达基因的蛋白质相互作用网络，并使用Cytoscape 软件对蛋白互作的结果进行可视化（图8）。随后使用软件中的MCODE 插件筛选子模块，并设置MCODE 评分＞10 为标准。经过筛选后，我们从中获得了32 个基因（平均MCODE 得分=10.129），包括32 节点（Node）及157 个节点间的连线（Edge），其中一个节点代表一个蛋白质，节点间的连线代表蛋白间的相互作用。然后使用软件中的CytoHubba 插件，通过12 种算法筛选连接度最高的前10 位的基因作为关键基因（图9），分别为APP、DYNC1H1、甲酰肽受体2（formyl peptide receptor 2，FPR2）、GNB1、毒蕈碱型胆碱受体M2（cholinergic receptor muscarinic 2，CHRM2）、CD59、ARPC2、GNG7、AP2M1、HSP90B1。

3 讨论

HCM 可导致心律失常性猝死、心力衰竭和心房纤颤（可致栓塞性卒中）等不良后果，是年轻人包括训练有素的运动员猝死的最常见原因之一［6］，其治疗方式现已有了较大的进步，但其发病机制仍未完全明确，许多人的治疗效果并不理想，预期寿命或实质性症状没有大幅度降低［7］。关于HCM的基因研究始于20多年前，现已发现上千种肌节相关蛋白基因如MYH7、MYPBC3 和TNNT2 等基因的突变与疾病的发生相关［8］。基因检测可在临床疾病发病前识别无症状的HCM 患者，从而对疾病的诊断、治疗和预后进行指导，因此，基因检测已成为欧洲和北美指南中I 类推荐的辅助检查［9-11］。由此可见，对HCM 的基因谱进行分析，从而对HCM 的发病机制进行分子层面的探索，这对疾病的预防、治疗及预后有着极其重要的作用。

本研究通过从GEO 数据库中下载GSE36961数据集并使用R 语言的“limma 程序包”进行分析后，共获得与HCM 相关的差异表达基因8 002 个。这些数据表明，HCM 患者的心肌组织与正常人的心肌组织间的基因表达确实存在差异。

为系统性展现并解释疾病的分子网络机制，本研究利用了R 语言中的“WGCNA 程序包”进行网络图谱构建，从中发现了与HCM 相关的青色模块与紫红色模块。通过对青色模块内的基因进行功能富集分析，本研究发现其生物过程、分子功能及信号通路主要富集在能量代谢中，细胞组分主要富集在线粒体中。既往研究发现，在心肌病理性肥大的发展过程中，心肌细胞重构了能量代谢的产生机制，减少了葡萄糖的氧化，增加了糖酵解和其他形式的代谢［12］，这种细胞代谢的改变直接诱导或促进了心肌病理性肥大的发生［13］。此外，供应心肌细胞能量的减少导致了心肌细胞的死亡和纤维化，从而导致心肌细胞从适应到心肌肥大甚至心力衰竭的转变［14］，而线粒体功能障碍是这一过程公认的潜在机制［15-16］。这些都和本研究所发现结果相符合。

通过对紫红色模块内的基因进行功能富集分析，本研究发现其生物过程主要富集在血管形成中。既往研究表明，毛细血管的数量会对心肌能量的供应造成影响，从而成为影响心肌病理性肥大的重要因素［17］。血管内皮生长因子（vascular endothelial growth factor，VEGF）是维持心肌新生毛细血管数量的重要血管生成分子，VEGF 的缺失会导致心肌血管生成的减少和心脏功能的受损［18］。目前VEGF 信号通路抑制剂已被临床上批准用于治疗各种类型的癌症，然而，这些抑制剂的使用可能导致心肌病等心血管相关毒性疾病的发生［19］。在心肌生理性肥大中，新生毛细血管的数量随着心肌细胞的生长而增加，从而为心肌提供足够的营养和氧气，但在心肌病理性肥大时，新生毛细血管数量和冠状动脉血流储备不足，不能支持心肌生长，从而导致心肌轻度缺氧和营养不足［20-21］。这些也与本研究所发现结果相一致。

本研究还发现，紫红色模块内基因的细胞组分、分子功能和信号通路主要富集在细胞外基质成分的形成当中。本研究还通过STRING 数据库筛选出前10 位连接度最高的基因，其中FPR2 被报道与心肌的炎症反应及免疫激活相关［22-23］，CHRM2 被报道与自身免疫性心肌病相关［24］。既往研究显示，心肌病理性肥大的发生与炎症信号因子的释放和免疫细胞的激活有关。而炎症信号因子在心肌细胞和非心肌细胞中所引发的反应各不相同，比如白细胞介素-1β、肿瘤坏死因子-α和白细胞介素-6 等炎症信号因子，在心脏的成纤维细胞中，它们所引发的反应主要表现为抑制心肌细胞增殖，减少基质合成，以及增加基质金属蛋白酶的活性；在心肌细胞中，则表现为诱导心肌细胞肥大，并引起心肌细胞的凋亡；而在免疫细胞中，则表现为促进炎症的发生，造成心肌细胞的损伤［25-26］。由此可见，在心脏损伤和肥厚性重构的过程中，心肌炎症的发生可以诱导心肌病理性肥大和心肌纤维化，最终导致心肌细胞死亡的不良结局［27］。以上亦与本研究所发现结果相符合。

此外，本研究发现年龄或性别与模块之间没有显著相关性。这与既往流行病学调查结果相类似：HCM 与性别、年龄无明显相关性，无论男女均可发生，亦可在任何年龄出现临床表现，儿童和成人均可能发生猝死［28-29］。

综上所述，本研究通过生物信息学算法再分析了GEO 数据库中106 例HCM 及39 例对照组人类心肌样本，一方面充分利用数据库资源，减少了取样和测序的成本，另一方面通过更新的系统生物信息学算法也可再挖掘更多的分子机制信息。本研究通过分析结果识别出与HCM 相关的两个模块，并筛选出相应的关键基因，印证了HCM 的发病机制与能量代谢、血管形成及炎症反应相关，为HCM 的进一步治疗提供思路，从而使患者能够最大程度从治疗中获益。同时，识别出的关键基因中，除FPR2 及CHRM2 外的其余基因与HCM 的关联虽尚未见报道，但其对HCM 发病的分子机制的深层次分析具有一定的价值，这为HCM 分子机制的深入研究提供了一定的理论依据，这些基因的作用仍需在未来的细胞及动物实验中得到进一步的验证。