Stanford A型主动脉夹层基因表达谱分析

高原,李召水,李磊,孙展发,池一凡,生伟

(青岛大学附属青岛市市立医院心脏外科,山东 青岛 266071)

Stanford A型主动脉夹层(ATAAD)为一种常见的严重威胁人类生命健康的主动脉灾难性事件，主要发生在70～80岁年龄人群[1]。在过去的40年中，随着诊疗技术及治疗手段的进步，ATAAD不再是一种死亡后诊断疾病[2]。美国一项研究结果显示，ATAAD发生的频率几乎为腹主动脉瘤破裂的3倍[3]。且三级护理中心记录的ATAAD研究结果表明， ATAAD仍为通过区域快速运输系统转移的最频繁的急诊疾病[4-5]。因此，针对ATAAD发生风险的早期筛查尤为重要。

通过生物信息学手段进行生物大数据分析已经在多种生物学领域被应用。由于实验条件的限制性、样本收集的可行性等，分析数据库中已有数据可为研究疾病发生的分子机制提供便捷途径[6-7]。研究显示，ATAAD的发生与基因突变及表达异常具有密切关系，各种结缔组织基因的突变使个体易发生主动脉夹层[8]。本研究通过分析GEO数据库中ATAAD病人基因表达芯片测序大数据，进一步深入研究ATAAD发病的基因水平调控机制。现将结果报告如下。

1 资料和方法

1.1 基因微阵列数据检索

ATAAD样本的基因芯片来自于GEO(Gene Expression Omnibus)数据库(http://www.ncbi.nlm.nih.gov/geo/)。将Aortic dissection输入到GEO数据库检索域，结果显示共有158项研究。去除77项小鼠样本研究，对剩余81项人类样本研究进行进一步筛选。由于本文目的为研究ATAAD病人基因表达整体模式，需要搜寻得到一个芯片或测序整体信息，因此根据数据库条目类型，选择了包含11项相关研究的数据集类型。随后进一步筛选得到其中3项(GSE147026、GSE52093、GSE98770)研究为ATAAD病人mRNA相关的测序分析，由于GSE147026样本标注不明确，本文最后选择符合研究标准的2项即GSE52093与GSE98770进行ATAAD基因表达谱分析。相关芯片信息见表1。

表1 ATAAD芯片相关资料

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

下载得到原始数据后，首先通过R语言的Bioconductor package (BiocGenerics, Biobase, Parallel)对其进行噪声去除、分位数归一化等。随后，使用R语言(Limma包)于芯片数据中筛选ATAAD和正常样本之间的DEGs。筛选DEGs的阈值设定为表达上调或下调2倍，且P值<0.05。

1.3 火山图对DEGs进行可视化分析

将所得到的每一样本中每一基因的代表表达水平的数值、每一基因在实验样本组中的相对表达值及P值，按照R语言中volcano要求制作对应矩阵列表，按照试验需求设置阈值，随后通过R语言进行DEGs可视化分析。

1.4 热图(heatmap)对DEGs进行可视化分析

首先将分析得到的DEGs、测序样本及每一样本中每一基因对应表达水平的数值，按照R语言中pheatmap包要求制作对应矩阵列表，随后通过R语言(pheatmap包)对基于mRNA表达水平的组样本进行可视化层次聚类分析。

1.5 GSE52093和GSE98770中co-DEGs分析

GSE52093和GSE98770中的co-DEGs通过Draw Venn Diagram线上数据库进行分析(http://bioinformatics.psb.ugent.be/webtools/Venn/)。将待分析基因列表上传到数据库，随后即显示韦恩图及相关的共同基因列表。

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

将基因列表上传至DAVID生物信息学资源6.8数据库(https://david.ncifcrf.gov/)，随后显示GO(biological process, cell component, molecular function, KEGG)结果，将其下载为txt文件。最后通过R语言(ggplot2包)可视化GO结果。

1.7 蛋白质间相互作用关系网络分析

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

2 结 果

2.1 ATAAD病人基因表达变化分析

根据筛选条件在GSE98770数据集中得到的DEGs中包含上调基因239个(0.720%)，下调基因575个(1.733%)，无显著性差异的基因32 254个(97.257%)(图1A、B)。在GSE52093数据集中以相同的条件进行分析，最终得到上调基因315个(0.654%)，下调基因308个(0.640%)，无显著性差异的基因47 384个(98.49%)(图1C、D)。以上结果表明，在ATAAD疾病发展的过程中，有1.3%～2.5%基因表达水平发生改变。

A：火山图显示GSE98770数据集中ATAAD病人与正常人相比具有表达差异的基因，其中红色表示ATAAD中表达上调的基因，绿色表示ATAAD中表达下调的基因，灰色表示为表达无显著性差异的基因(横坐标表示基因表达倍数的对数值，纵坐标表示P值对数值)；B：热图显示GSE98770数据集中被测样本DEGs的表达趋势，其中蓝色代表低表达，红色代表高表达：C：火山图显示GSE52093数据集中ATAAD病人与正常人相比具有表达差异的基因，其中红色表示为ATAAD中表达上调的基因，绿色表示为ATAAD中表达下调的基因，灰色表示为表达无显著性差异的基因(横坐标表示基因表达倍数的对数值，纵坐标表示P值对数值)；D：热图显示GSE52093数据集中被测样本DEGs的表达趋势，其中蓝色代表低表达，红色代表高表达。

随后对 GSE98770与GSE52093数据集中共同的DEGs进行分析，结果显示，共同上调2倍基因27个(图2A)，共同下调2倍基因67个(图2B)，为方便描述，后文将其命名为co-DEGs。

左图：GSE52093与GSE98770数据集中共同表达上调基因的数目(27个)；右图：GSE52093与GSE98770数据集中共同表达下调基因的数目(67个)。

2.2 ATAAD病人中DEGs生物学功能

应用GO分析对前面得到的94个co-DEGs的生物学功能(biological process)、细胞内成分(cellular components)、分子功能(molecular function)及KEGG进行分析，结果显示，ATAAD病人样本中DEGs主要涉及的生物学过程为调节血管收缩、细胞氯离子稳态、心室心肌细胞动作电位、连环蛋白导入细胞核的正向调节、肌动蛋白丝聚合的负性调节、钠离子跨膜转运体活性的调控及平滑肌收缩的正向调节等(图3A)。DEGs涉及的分子功能为钙诱导的钙释放活性、ryanodine敏感的钙释放通道活性、肌酸激酶活性、泛素蛋白连接酶活性、钙离子结合及钠通道调节器活性等(图3B)。DEGs所处的细胞成分为肌纤维膜、肌动蛋白骨架、细胞质、肌肉神经接点、光滑型内质网、黏着斑及肌浆网膜等(图3C)。DEGs参与的KEGG包括血管平滑肌收缩、精氨酸和脯氨酸代谢、心肌细胞的肾上腺素能信号、缩宫素信号通路、cGMP-PKG信号通路、心律失常性右心室心肌病、PPAR信号通路、肥厚型心肌病、钙信号通路以及扩张型心肌病等通路(图3D)。GO分析结果表明，ATAAD疾病的发生可能与心肌细胞的膜内外转运调控、平滑肌收缩异常调控以及心肌细胞疾病有直接的关系。

A：co-DEGs参与的细胞内生物学过程；B：co-DEGs在细胞内的分子功能；C：co-DEGs涉及的细胞组分；D：co-DEGs参与的通路。纵坐标表示分析得到的生物学过程、分子功能、细胞组分或通路，横坐标表示P值的对数值。

2.3 蛋白质互作网络分析ATAAD中DEGs相互关系

对前面94个co-DEGs的蛋白质相互作用关系进行分析，结果显示，互作网络为以RYANODINE受体2(RYR2)、心肌素(MYOCD)及平滑肌动蛋白2(ACTG2)等为核心的关系网络，以及以血管紧张素Ⅱ受体Ⅰ型(AGTR1)和ASB2为核心的关系网络(图4)。此外，两数据集中均分析得到共同上调较为显著的类血管生成素4(ANGPTL4)、G蛋白信号转导调节因子2(RGS2)以及溶质载体家族16成员3(SLC16A3)等。

图4 ATAAD中DEGs之间相互关系的蛋白质互作关系

对蛋白互作网络中心以及表达水平较为突出的基因ACTG2、AGTR1、RYR2、MYOCD、RGS2、ANGPTL4、SLC16A3、精胺氧化酶(SMOX)表达水平进行分析，结果显示，GSE98770及GSE52093数据集中的ATAAD病人样本中ACTG2、AGTR1、RYR2、MYOCD相对表达水平均低至0.4(图5A)，ANGPTL4、RGS2、SLC16A3、SMOX相对表达水平均最高可至12倍(图5B)，表明以上基因的异常表达与ATAAD疾病的发生具有重要关系。此外，关键基因所涉及的KEGG多为平滑肌收缩调控、肥厚型心肌病、心肌细胞的肾上腺素能信号、缩宫素信号通路、钙离子代谢等(图5C)，提示ATAAD疾病的发生归因于平滑肌收缩异常调控及心肌细胞膜运输异常调控。

图5 关键差异基因相对表达水平及其GO分析

3 讨 论

本文研究通过生物信息学的方法，对GEO数据库GSE98770和GSE52093数据集中ATAAD病人体内基因表达数据进行分析，结果显示ATAAD病人中1.3%～2.5%的基因表达水平发生了改变。

GSE98770与GSE52093虽为ATAAD病人DEGs测序分析，但二者的不同之处在于GSE98770对病人升主动脉夹层中的内膜层组织进行测序分析，而GSE52093对病人升主动脉夹层整体样本进行测序分析。由于两个数据集中病人来源于不同国家地区，且测序平台不同，取样部位不同，因此通过对两个数据集中分析得到的DEGs进行交集分析后共筛选得到94个co-DEGs，表明ATAAD的发生与基因表达模式的改变具有重大关联。

GO分析通常用于研究DEGs生物学功能的富集情况，以研究疾病发生过程中所发生改变的生物学过程。本文GO结果显示，co-DEGs生物学功能显著富集于心肌细胞的膜内外转运调控、平滑肌收缩异常调控的生物学过程。本文蛋白互作关系分析筛选得到位于网络核心的关键基因ANGPTL4、AGTR1、SLC16A3、MYOCD、ACTG2、RGS2、SMOX、以及RYR2，且这8个关键基因在ATAAD病人样本中表达水平均发生改变，表明ACTG2、AGTR1、RYR2、MYOCD、ANGPTL4、RGS2、SLC16A3以及SMOX可能对于ATAAD的发生发展过程具有关键作用。ATAAD的发生与这些关键基因表达水平的改变或呈一定相关性。

已有研究显示，ACTG2、AGTR1、ANGPTL4、MYOCD、RYR2、RGS2、SLC16A3以及SMOX均参与人类多种疾病的调控。ACTG2参与平滑肌细胞的调控[9]，且该基因突变与多种肠道疾病的发生具有重要关系[10-12]。此外，AGTR1编码1型受体，被认为介导了血管紧张素Ⅱ的主要心血管作用，并可能在缺血或梗死心肌血流恢复后再灌注性心律失常中发挥作用，以及介导低氧对早期胎盘血管生成的影响[13]。相关研究结果显示，RYR2的突变与应力诱导的多型室性心动过速和心律失常性右心室发育不良有关，且该基因的异常表达参与儿茶酚胺敏感性室速的诊断与治疗[14-15]。MYOCD参与人胸膜间皮细胞间充质转化、心肌细胞死亡及分化[16-17]，基因ANGPTL4突变与代谢综合征及脑动静脉畸形风险显著相关[18-19]。RGS2参与多种人类神经性疾病的发生[20-22]。SLC16A3在细胞代谢和组织传递中发挥重要作用，并参与人类多种癌症的调控[23-24]。SMOX在癌症的发生中具有一定作用[25-26]。

本研究利用生物信息学技术及其方法探讨了ATAAD病人样本DEGs，分析可能与ATAAD疾病发生相关的DEGs以及相关生物学途径。结果显示,ATAAD的发生与生物体内多种基因表达紊乱并造成相应调控途径改变相关，研究DEGs可从基因水平解释ATAAD的发病机制，并为ATAAD的临床治疗及诊断提供新的思路。后续的研究中将进行DEGs体外功能验证以进一步研究相关基因与ATAAD发病机制的关系。