1. 引言
冠心病、缺血性中风和外周动脉疾病等多种心脑血管疾病是目前全球死亡率和致残疾率较高的疾病[1]。动脉粥样硬化(atherosclerosis, AS)是这些疾病的病理基础。AS是一种发生于血管壁的慢性炎症性病理变化，以脂质沉积和免疫细胞浸润为特征[2]。衰老是生命过程中机体发生的随时间推移而形成的退行性变化。机体衰老伴随着生理过程和解剖结构复杂性的丧失[3]。细胞衰老是指细胞在一系列内源性和外源性刺激源的作用下经历独特的表型改变，导致细胞分裂停止，细胞周期通常被永久阻滞在G0或G1期[4]。在机体衰老过程中，衰老细胞在心血管组织和器官上逐渐累积，机体发生组织修复能力的丧失和功能障碍[5]。因此，减少细胞衰老，进而延缓衰老相关心血管疾病的发生具有重要意义。随着研究的深入，研究者发现血管中的多种衰老细胞与AS病理生理学变化密切相关，同时衰老细胞所分泌的分泌表型(senescence-associated secretory phenotype, SASP)也会导致AS斑块的进展和失衡[6]。因此，细胞衰老逐渐被作为动脉粥样硬化性心血管疾病的主要危险因素之一。
本研究运用生物信息学方法，对源自基因表达谱数据库(GEO)的AS相关数据集进行差异基因表达分析，筛选AS形成的差异表达基因(differentially expressed genes, DEGs)，联合衰老基因和免疫浸润分析，为靶向衰老细胞延缓衰老相关心血管疾病的发生提供理论依据。
2. 资料与方法
2.1. 数据资源
本研究的三组转录组数据来自于基因表达数据库(GEO) (http://www.ncbi.nlm.nih.gov/geo)。其中GSE21545芯片数据包括223例接受动脉内膜切除术的患者的外周血单核细胞基因表达量和颈动脉斑块基因表达量，为方便研究，我们选取了其中126例颈动脉斑块基因表达量。GSE43292芯片数据包括32例高血压患者的颈动脉斑块基因表达量和配对的临近斑块正常颈动脉组织的基因表达量，两组数据集分别采用Affymetrix Human Gene 1.0 ST Array和Affymetrix Human Genome U133 Plus 2.0 Array平台。GSE100927芯片数据集作为外部验集。307个人类衰老基因(aging-related genes, ARGs)来自人类衰老基因组资源库(https://genomics.senescence.info/)。
2.2. 数据处理
将GSE21545和GSE43292中的数据合并，并使用R包“limma”和“sva”去除批次效应，减少数据集之间由于实验误差导致的差异，让多个数据集重新组合在一起，使得下游分析可以只考虑生物学差异因素。
2.3. DEGs的筛选及可视化
应用R语言的limma包对合并的数据集中AS组和对照组的基因表达量进行分析，应用倍比法(fold-change, FC)和假发现率(false-discovery rate, FDR)筛选GEGs，为扩大筛选条件：仅限PDR值 < 0.05。对GEGs进行聚类热图及火山图展示[7]。
2.4. DEGs的功能分析及通路分析
应用R语言的clusterProfiler包对GEGs进行GO基因功能分析及KEGG的Pathway富集分析，并进行可视化展示。
2.5. 筛选差异表达的衰老相关基因(Differentially Expressed Aging-Related Genes, DE-ARGs)，进行富集分析及蛋白相互作用网络
将DEGs与ARGs取交集，得到与AS相关De-ARGs，进行GO基因功能分析及KEGG的Pathway富集分析。使用STRING (http://string-db.org)进行蛋白相互作用网络分析，构建De-ARGs的蛋白相互作用网络，选择综合分数大于0.4的主要蛋白相互作用网络(protein interaction network, PPI)，将得出的PPI数据导入Cytoscape软件进行进一步分析。运用Cytoscape软件的cytoHubba插件，分别采用MCC、MNC、EPC、Degree、Closeness、Radiality六种方法，选取各自前10位的Node基因取交集作为hub基因。
2.6. 免疫细胞浸润程度的评估
免疫细胞浸润的分析在了解疾病进展和治疗反应方面具有重要的驱动作用。利用23个免疫基因集构建了GSEA基因集的单样本基因集富集分析(ssGSEA)。使用R语言GSVA包的ssGSEA函数评估免疫细胞浸润比例，比较AS样本与对照组样本免疫细胞浸润差异。此外，我们检测了每个hub基因与浸润的免疫细胞的关系，并绘制热图进行可视化展示。
2.7. 基因集变异分析(Gene Set Variation Analysis, GSVA)
为进一步了解hub基因在AS组中的整体趋势，我们对hub基因在AS组的表达进行GSVA分析，采用http://software.broadinstitute.org/gsea/msigdb数据库中的C2数据集作为背景基因。
2.8. ceRNA网络的构建
利用miRanda、miRDB、TargetScan三个数据库预测hub基因的miRNA，同时应用spongeScan数据库预测miRNA的lncRNA。将结果导入Cytoscape软件构建hub基因调控的ceRNA网络。
2.9. 外部数据集验证
利用GSE100927芯片数据进一步验证hub基因在两组中的表达，并绘制箱线图进行可视化展示。
2.10. 统计学方法
数据统计分析由R4.3.0完成。使用Wilcox检验，评估两组之间的统计学差异。Spearman相关性用于分析hub基因的表达水平与免疫细胞之间的关系。P < 0.05具有统计学意义。
3. 结果
3.1. 数据处理
选取GSE21545芯片数据样本中126例颈动脉斑块基因表达量，GSE43292芯片数据样本中32例高血压患者的颈动脉斑块基因表达量和配对的临近斑块正常颈动脉组织的基因表达量，根据平台注释信息，进行探针–基因名转换。根据基因名将两组数据合并，合并后的数据集包含158例AS样本和32例对照组样本。使用R包“limma”和“sva”去除批次效应，减少数据集之间由于实验误差导致的差异，基于PCA，如图1(A)所示，两个数据集合并后各自独立，经去批次效应处理后，两个数据集合并后具有较好融合(图1(B))。应用limma包的normalize Between Arrays函数，对表达矩阵进行归一化处理，消除奇异样本数据导致的不良影响。
3.2. DEGs的筛选及可视化
应用R语言的limma包对标准化处理后的合并数据集的158个AS化样本和32个对照组样本的基因表达量进行分析，筛选PDR值 < 0.05的差异表达基因3817个，其中上调基因1821个，下调基因1996个，对GEGs进行火山图展示(图2(A))，红色为上调基因，蓝色为下降基因。对DEGs绘制聚类热图(图2(B))，在图中可看出，合并数据集中动脉粥样硬化组与对照组在基因表达上具有较好的区分。
3.3. DEGs的功能富集分析
对DEGs的GO富集分析表明(图3(A))，DEGs的生物学过程(biological process, BP)主要富集在细胞因子产生的正调控，对外部刺激反应的正调节，细胞粘附的正调控；分子生物学功能(molecular function, MF)主要富集在肌动蛋白结合、GTP酶调节活性、核苷–三磷酸酶调节活性；细胞学组分(cellular components, CC)主要富集在细胞–基底连接、局灶性粘附、细胞前缘。KEGG通路分析显示，脂质与动脉粥样硬化、趋化因子信号通路、细胞粘附分子被显著富集(图3(B))。
3.4. 筛选De-ARGs，进行富集分析
将DEGs与ARGs取交集，得到与AS相关差异表达ARGs (De-ARGs)，如图4(A)，共筛选出88个De-ARGs。其中46个上调基因，42个下调基因。对88个De-ARGs进行GO富集分析表明(图4(B))，BP主要富集于凋亡信号通路，的调控，对肽的细胞反应，对氧化应激的反应等。MF主要富集于蛋白丝氨酸/苏氨酸/酪氨酸激酶活性，蛋白磷酸酶结合，磷酸酶结合等。CC主要富集于囊泡腔，膜微区，PML小体。KEGG通路分析显示，PI3K-Akt信号通路，FoxO信号通路，JAK-STAT信号通路等被显著富集(图4(C))。
(A) (B)
Figure 1. Data processing
图1. 数据处理
(A) (B)
Figure 2. Screening of DEGs. (A) Merge volcano plots of differentially expressed genes in the dataset; (B) Heat map of differential expression of top 20 genes between atherosclerosis group and control group in data set
图2. DEGs的筛选。(A) 合并数据集中差异表达基因的火山图；(B) 数据集中动脉粥样硬化组与对照组差异表达前20位基因的热图
(A)
(B)
Figure 3. Functional enrichment analysis of DGEs. (A) GO enrichment analysis bubble diagram of DGEs; (B) KEGG enrichment analysis Circos diagram of DEGs
图3. DGEs的功能富集分析。(A) DGEs的GO富集分析气泡图；(B) DEGs的KEGG富集分析Circos图
(A) (B)
(C)
Figure 4. Screening and enrichment analysis of De ARGs. (A) Veen plot of intersection genes between DEGs and ARGs; (B) GO enrichment analysis of De ARGs; (C) KEGG enrichment analysis of De ARGs
图4. De-ARGs的筛选及富集分析。(A) EGs与ARGs交集基因的Veen图；(B) De-ARGs的GO富集分析；(C) De-ARGs的KEGG富集分析
3.5. 蛋白质相互作用和核心基因分析
我们将87个De-ARGs输入在线STRING数据库，构建了De-ARGs的PPI网络(图5(A))，使用cytoHubba插件，分别采用MCC、MNC、EPC、Degree、Closeness、Radiality六种方法选择前10位的Node基因取交集作为hub基因(图5(B))。最终IGF1，GRB2，MYC，EGFR，PTEN5个基因被鉴定为hub基因。其中，IGF1，GRB2，MYC为上调基因，EGFR，PTEN为下调基因。
(A)
(B)
Figure 5. PPI network interactions and hub gene analysis. (A) The PPI network of DE ARGs; (B) Six methods are used to calculate the Upset plot of the intersection of the top 10 Node genes
图5. PPI网络相互作用和hub基因分析。(A) DE-ARGs的PPI网络；(B) 六种方法计算前10位Node基因取交集的Upset图
3.6. 免疫细胞浸润程度的评估
我们使用R语言GSVA包的ssGSEA函数评估免疫细胞浸润比例，比较AS样本与对照组样本免疫细胞浸润差异。结果显示(图6(A))，与对照组相比，AS组样本除CD 56 dim natural killer cell和Type 2 T helper cell外21种免疫细胞均升高，差异具有统计学意义(P < 0.05)。同时我们比较了5个hub基因与免疫细胞浸润的相关性(图6(B)~(F))，分析显示IGF1，MYC与所有免疫细胞成正相关，与之相反，EGFR与所有免疫细胞均成负相关。而GRB2与Type 2 T helper cell成负相关,与其他免疫细胞均成正相关。PTEN与绝大多数免疫细胞成负相关，Type 2 T helper cell除外。
(A)
(B) (C)
(D) (E)
(F)
Figure 6. Evaluation of immune cell infiltration degree. (A) The difference in immune cell infiltration between AS samples and control group samples; (B)~(F) The correlation between hub genes and immune cell infiltration
图6. 免疫细胞浸润程度评估。(A) AS样本与对照组样本免疫细胞浸润差异；(B)~(F) hub基因与免疫细胞浸润的相关性
3.7. hub基因的GSVA富集分析
GSVA富集分析是一种特殊类型的基因集富集方法，通过将基因集在不同样品间的表达量矩阵转化成基因集在样品间的表达量矩阵，从而在生物信息学角度解释导致表型差异的原因。通过GSVA分析(图7)，我们发现IGF1、GRB2和MYC在各自高表达组中，均出现自然杀伤细胞介导的细胞毒性、Toll样受体信号通路、NOD-like受体信号通路等多个通路均上调。而基底细胞癌通路、牛磺酸和低牛磺酸代谢途径和丙酸代谢通路均下调。与之相反，EGFR、PTEN作为负向调节基因，在均在丙酸代谢通路上调，EGFR在而基底细胞癌通路、牛磺酸和低牛磺酸代谢途径也明显上调。可见枢纽基因间存在丰富交叉富集通路。
(A)
(B)
(C)
(D)
(E)
Figure 7. GSVA enrichment analysis of hub genes
图7. hub基因的GSVA富集分析
3.8. ceRNA网络的构建
我们利用miRanda、miRDB、TargetScan三个数据库预测hub基因的上游miRNA，然后应用spongeScan数据库预测miRNA的lncRNA。将结果导入Cytoscape软件构建hub基因调控的ceRNA网络。如图8所示，mRNA与miRNA、lncRNA之间具有一对多和多对一的对应关系，可见它们之间可能具有较复杂的调控关系。
3.9. 外部数据集验证
在GSE100927数据集中，我们验证了hub基因的表达。如图9(A)~(C)所示，与训练数据集相同，与对照组相比，AS组样本IGF1，GRB2，MYC表达显著升高，PTEN则显著下调(图9(D))，差异具有统计学意义(P < 0.05)。而与训练集相反，GSE100927数据集中EGFR在AS组表达高于对照组(图9(E))。
4. 讨论
动脉粥样硬化(atherosclerosis, AS)是一种以脂质代谢紊乱为基础的病理改变。病变动脉内膜可见脂质沉积、纤维组织增生和钙化，随后动脉壁中层逐渐出现变性和钙化，进而出现动脉壁增厚，弹性降低，血管变窄，常导致斑块破裂、出血、血栓形成等并发症，是冠心病、脑梗死、外周血管病的主要原因[8]。
我国已进入老龄化社会。据估计，到2030年，65岁以上的人口将占总人口的20% [9]。导致该年龄组人口死亡的首要因素是心血管疾病，占比高达40%，针对心血管疾病治疗的费用也将增加两倍[10] [11]。在机体衰老过程中，衰老细胞在心血管组织和器官上逐渐累积，机体发生组织修复能力的丧失和功能障碍。随着研究的深入，研究者发现血管中的多种衰老细胞与AS病理生理学变化密切相关，同时衰老细胞所分泌的SASP也会导致动脉粥样硬化斑块的进展和失衡[6]。因此，细胞衰老逐渐被作为动脉粥样硬化性心血管疾病的主要危险因素之一。
Figure 8. Construction of ceRNA network
图8. ceRNA网络构建
(A) (B)
(C) (D)
(E)
Figure 9. Expression of hub gene in GSE100927 dataset
图9. hub基因在GSE100927数据集中的表达
近年来研究表明，靶向衰老细胞有望延缓衰老相关心血管疾病的发生[5]，为延缓衰老和防治心血管疾病提供新机遇。本研究运用生物信息学方法，对源自基因表达谱数据库(GEO)的AS相关数据集进行差异基因表达分析，筛选AS形成的差异表达基因(differentially expressed genes, DEGs)，联合衰老基因和免疫浸润分析，为靶向衰老细胞延缓衰老相关心血管疾病的发生提供理论依据。经本研究发现，IGF1，GRB2，MYC，EGFR，PTEN5个基因为与衰老相关的AS关键基因，其中，IGF1、GRB2、MYC为上调基因，EGFR，PTEN为下调基因。经外部数据集验证，发现IGF1，GRB2，MYC，PTEN的变化趋势与训练集相同，在衰老相关AS中发挥重要作用。
PTEN是人第10号染色体缺失的磷酸酶及张力蛋白同系物，不仅调控细胞的生长和存活，而且作为一种多功能蛋白，在很多细胞进程中发挥广泛的作用。PTEN通过其磷酸酶活性，负向调控PI3K/AKT信号通路，进而控制细胞的生长、增殖和凋亡。PTEN去磷酸化PIP3 (磷脂酰肌醇三磷酸)，抑制AKT的激活，从而阻止下游信号传导，这对于维持细胞正常功能和防止过度增殖至关重要[12] [13]。PTEN还能够影响细胞周期的进程，其通过调节Cyclin D1和p27等细胞周期蛋白，控制细胞从G1期向S期的过渡，防止细胞异常增殖[14]。本研究显示，与对照组相比，AS组PTEN表达明显下降。免疫浸润分析显示，PTEN与除Type 2 T helper cell之外的其他免疫细胞均成负相关。PTEN作为1997年发现的第一个具有双重磷酸酶活性的肿瘤抑制基因，同时具有脂质磷酸酶和蛋白磷酸酶双重活性，能够调控细胞的生长、凋亡等过程，并且在癌症发生、心肌缺血再灌注、心肌肥厚等疾病发生过程中也具有调控作用，对于血管内皮细胞生物学行为、心脏发育、胰岛素信号转导等具有重要意义[15] [16]。AS的发生与细胞免疫炎症、血管内皮细胞凋亡等相关[17] [18]。PTEN具有抗炎作用，能够通过调控中性粒细胞的趋化减弱组织炎症，并且能够抑制血管平滑肌细胞的迁移，延缓动脉粥样硬化的发生[19] [20]。Tsoyi [21]等的研究发现，PTEN可以通过PI3K/Akt/GSK-3/GATA-6信号通路选择性抑制VCAM-1的表达。本研究通过GSVA分析显示，在AS组中，PTEN低表达组B细胞受体信号通路，T细胞受体信号通路，天然杀伤细胞介导的细胞毒性，细胞凋亡等通路明显上调。说明PTEN下调，通过激活多通路参与动脉粥样硬化的发生发展。
MYC基因是较早发现的一组癌基因，包括C-myc，N-myc，L-myc，分别定位于8号染色体，2号染色体和1号染色体。本研究显示，AS组MYC较对照组明显上调。C-myc基因受血小板源性生长因子(platelet derived growth factor, PDGF)等多种有丝分裂原刺激后呈快速诱导方式表达，参与细胞增殖调控，是细胞生长过程中的重要调节基因。C-myc的产物是核内蛋白质，与细胞G0-G1期有关，参与DNA复制内膜损伤、高脂血症等多种因素，可引起血管壁PDGF基因的异常表达，而诱导C-myc基因表达核内蛋白质，与DNA结合，条件细胞内其他基因表达，产生和分泌多种生长因子，再通过内分泌和旁分泌作用，使细胞增生更明显[22]。本研究通过芯片GSE100927数据也证实，MYC在AS组表达上调，参与动脉粥样硬化的形成。免疫浸润分析显示，MYC与所有免疫细胞成正相关，可能通过免疫调控参与动脉粥样硬化。同时，GSVA显示，MYC高表达组，Toll样信号通路，NOD-like受体信号通路，细胞凋亡通路，T细胞受体信号通路，B细胞受体信号通路等参与动脉粥样硬化的通路均有显著富集。
GRB2即生长因子受体结合蛋白2 (growth factor receptor-bound protein 2, GRB2)，是一种桥接蛋白，通常在许多类型的细胞中正常表达。GRB2调节的信号通路对于维持细胞增殖、生长和分化等多种基本生理功能至关重要，其异常激活和异常表达与癌症细胞的增殖和侵袭以及各种其他慢性疾病的发生和发展有关[23] [24]。Brandon等人证实GRB2是AS病变形成和巨噬细胞摄取ox-LDL所必需的[25]。既往细胞和动物实验已经发现GRB2的表达水平与脂质的积累和炎症浸润的程度有关，敲低GRB2具有很强的抗AS作用[26]。本研究也证实，与对照组相比，AS组织中GRB2表达显著上调，同时，GRB2与除Th2辅助性T细胞外的其他免疫细胞成正相关。
胰岛素样生长因子-1 (insulin-like growth factor-1, IGF1)是胰岛素家族的一种单链多肽类生长因子，其结构与胰岛素具有高度同源性，具有细胞分化、增值功能和胰岛素样代谢作用。IGF1是目前已知的最有效的抗凋亡、促存活、抗氧化和细胞保护因子之一[27]-[30]。与既往研究不同，本研究显示，AS组IGF1明显上调，并且验证集数据也显示相同结果。产生该结果的具体原因尚不明确，我们准备下一步选取相关临床数据进一步分析。
本研究通过免疫分析显示，与对照组相比，AS组样本免所有疫细胞均升高，但CD56+自然杀伤细胞、Th2辅助T细胞差异无统计学意义。不同的免疫细胞在AS中既有抗炎作用，亦有促炎作用，且其参与炎症反应的作用机制不尽相同。与此同时，免疫细胞在AS的不同时期对斑块和炎症的促进和抑制作用不同。调节免疫细胞可能在未来抗AS及稳定易损斑块的治疗中具有巨大的潜力，但仍需要进一步的研究来充分阐明这些细胞在AS中的特定功能及其作用机制。
EGFR即表皮生长因子受体，是一种细胞膜结合受体，具有酪氨酸激酶活性，参与了主要信号通路的控制，包括细胞存活、增殖和迁移。EGFR过表达、自分泌配体刺激或构成型活性受体突变体[31] [32]可导致这种微调信号系统的失调，导致多种病理生理紊乱，促进癌症的发展。有研究显示，EGFR在高脂饲养的兔的外周血白细胞和动脉粥样硬化样斑块的巨噬细胞存在[33]，EGFR介导了体内对单核细胞的趋化反应和对单核细胞衍生的巨噬细胞的促有丝分裂反应。有报道显示，EGFR的参与对巨噬细胞的促动脉粥样硬化活性很重要，应用其抑制剂减少了促炎细胞因子的产生、脂质摄取和氧化应激[34]。本研究分析也显示，EGFR基因在AS组织中显著上调，然而在验证集GSE100927中，EGFR在AS组织中的表达则下调，需要下一步在基础实验中进一步挖掘EGFR参与动脉粥样硬化的证据。
本研究对分析计算得出的hub基因进行了GSVA分析，分析发现，自然杀伤细胞介导的细胞毒性、Toll样受体信号通路、NOD样受体信号通路等多个通路均上调。这些通路也被既往研究证实为参与动脉粥样硬化的重要通路。同时，不同hub基因间存在丰富交叉富集通路，说明hub基因通过多条通路参与AS的形成和发展。与此同时，我们利用miRanda、miRDB、TargetScan三个数据库预测hub基因的上游miRNA，然后应用spongeScan数据库预测miRNA的lncRNA，并构建了hub基因调控的ceRNA网络，它们之间具有一对多和多对一的复杂的调控关系。
本研究运用生物信息学方法，对源自基因表达谱数据库(GEO)的AS相关数据集进行差异基因表达分析，筛选AS形成的差异表达基因(differentially expressed genes, DEGs)，联合衰老基因和免疫浸润分析，为靶向衰老细胞延缓衰老相关心血管疾病的发生提供理论依据。
NOTES
*通讯作者。