基于加权基因共表达网络分析的动脉粥样硬化中铁死亡核心基因及其与免疫浸润细胞的关系研究

方柔柔，杨启帆，韩若冰，邬东东，孙娜，李娟，徐守竹，赵晶

作者单位：1.712046 陕西省咸阳市，陕西中医药大学公共卫生学院 2.712046 陕西省咸阳市，陕西中医药大学陕西省针药结合重点实验室

动脉粥样硬化（atherosclerosis，AS）是一种危害人类健康的慢性心血管疾病，每年约有2 000万人死于AS引起的疾病，且AS发病年龄呈年轻化趋势［1］。AS的主要病变是动脉壁存在脂质沉积并伴随平滑肌细胞增殖和纤维基质增生，逐渐形成动脉粥样硬化斑块［2］，从而引起冠状动脉狭窄和易损斑块形成，进而导致心血管事件风险增加［3］。在AS的发生发展过程中存在各种形式的细胞死亡，其中铁死亡与多种细胞病理过程相关［4］。铁死亡是一种新型非凋亡形式的调控性细胞死亡，其主要机制是氧化应激状态下铁代谢紊乱导致多种抗氧化系统异常，从而诱导脂质过氧化物在细胞内积聚［5］。研究显示，铁死亡在肿瘤、肝脏疾病、神经退行性疾病、内分泌系统疾病、心血管系统疾病中均是重要的治疗靶点［6-10］。因此，在AS发病率不断上升的背景下，研究AS中铁死亡的机制，寻找新的治疗靶点具有重要意义。

加权基因共表达网络分析（weighted gene coexpression network analysis，WGCNA）是一种将基因数据转换为共表达模块，从而探索基因与给定特征之间相关性的生物信息学方法。近年来，WGCNA已被广泛应用于鉴定各种疾病的生物标志物，并展现出良好的潜力［11-12］。本研究基于WGCNA筛选AS中铁死亡核心基因，并分析其与免疫浸润细胞的关系，以期为揭示AS与铁死亡之间的关系提供数据支撑，为探索新的AS治疗靶点提供依据。

1 资料与方法

1.1 数据来源

从基因表达综合（Gene Expression Omnibus，GEO）数据库（http://www.ncbi.nlm.nih.gov/geo/）下载AS转录组数据集GSE100927。该数据集是基于GPL17077平台（Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray 039381）的转录测序数据，共有104个样本。筛选出颈动脉样本41个，包含12个健康对照者的颈动脉样本（对照组）和29个AS患者的颈动脉样本（AS组）。从FerrDb数据库（http://www.zhounan.org/ferrdb/legacy/index.html）下载铁死亡相关基因（ferroptosis related genes，FRG），去除重复基因后，最终纳入564个FRG。

1.2 研究方法

1.2.1 数据预处理

利用R软件（版本4.2.3）对数据进行预处理，包括消除背景、标准化处理数据、构建表达矩阵等，同时去除缺失、重复及低表达的基因探针。

1.2.2 差异表达基因（differentially expressed genes，DEG）的筛选

使用R软件的limma包对GSE100927的41个颈动脉样本数据进行分析，主要筛选DEG，筛选标准为校正后P＜0.05和|log2FC|＞1。使用ggplot2包绘制火山图。

1.2.3 基因本体论（gene ontology，GO）功能富集分析和京都基因与基因组百科全书（kyoto encyclopedia of genes and genomes，KEGG）通路富集分析

使用R软件clusterProfiler包和org Hs.eg.db包对筛选出的DEG进行GO功能富集分析和KEGG通路富集分析，以P＜0.05为截断标准。

1.2.4 WGCNA

采用R包WGCNA对GSE100927的41个颈动脉样本的测序数据构建共表达网络，以获得与AS相关的基因模块。使用网络拓扑分析确定软阈值，并依据软阈值进一步计算拓扑重叠矩阵（topological overlap matrix，TOM），识别不同基因模块，计算每个基因模块与AS的相关系数及P值，筛选出与AS相关性最高的基因模块。取DEG、WGCNA中与AS相关性最高的基因模块和FRG的交集基因，并绘制其韦恩图。

1.2.5 蛋白质-蛋白质相互作用网络（protein-protein interaction networks，PPI）分析、核心基因的确定

通过在线String数据库（https://string-db.org/cgi/input.pl）构建交集基因的PPI网络图，设置最低互动分数为0.15，去除没有交互作用的基因后，通过Cytoscape软件（3.8.2版本）对网络进行可视化分析，利用cytoHubba插件的Degree算法计算出该网络中每个节点的Degree值，取Degree值排名前5的基因作为核心基因。

1.2.6 核心基因表达水平及验证

比较GSE100927中对照组和AS组核心基因表达水平。使用R软件中ggplot2包分析数据集GSE20129〔包括79个健康对照者的外周血样本（对照组）和56个AS患者的外周血样本（AS组）〕和GSE226790〔包括3个健康对照者的动脉样本（对照组）和3个AS患者的血管样本（AS组）〕中核心基因表达水平。

1.2.7 单样本基因集富集分析（single sample gene set enrichment analysis，ssGSEA）

从分子特征数据库（Molecular Signature Database，MSigDB）中下载“c2.cp.kegg.symbols.gmt”的注释基因集，进行ssGSEA，以进一步探究核心基因在AS进程中的调控作用。

1.2.8 免疫浸润分析

从CIBERSORT网站（http://cibersort.stanford.edu/）获取LM22文件，采用CIBERSORT算法估计22种免疫浸润细胞（初始CD4T淋巴细胞、静息CD4记忆T淋巴细胞、活化CD4记忆T淋巴细胞、静息肥大细胞、活化肥大细胞、静息树突状细胞、活化树突状细胞、静息自然杀伤细胞、活化自然杀伤细胞、中性粒细胞、嗜酸粒细胞、单核细胞、浆细胞、CD8T淋巴细胞、调节性T淋巴细胞、滤泡辅助性T淋巴细胞、γ-δT细胞、初始B细胞、记忆B细胞、M0型巨噬细胞、M1型巨噬细胞、M2型巨噬细胞）表达水平。应用ggcorrplot包分析核心基因表达水平与免疫浸润细胞表达水平的相关性。

1.3 统计学方法

利用R软件（版本4.2.3）进行数据分析。计量资料符合正态分布以（±s）表示，两组间比较采用独立样本t检验；计量资料不符合正态分布以M（P25，P75）表示，两组间比较采用Mann-WhitneyU检验。以P＜0.05为差异有统计学意义。

2 结果

2.1 DEG筛选结果

共筛选出了508个DEG，其中上调DEG 308个、下调DEG 200个，见图1。

图1 DEG分析的火山图Figure 1 Volcano map of DEG analysis

2.2 GO功能富集分析和KEGG通路富集分析结果

GO功能富集分析结果显示，DEG主要涉及的生物过程为参与细胞因子产生的正调节、细胞活化的正调控、白细胞介导免疫、外部刺激反应的正调节、骨髓白细胞活化、白细胞迁移、吞噬作用、参与免疫反应的白细胞活化、参与免疫反应的髓系细胞活化、吞噬作用的调节，主要涉及的细胞组分为分泌颗粒膜、含胶原蛋白的细胞外间质、内吞囊泡、膜筏、膜微区、分泌颗粒腔、三级颗粒、特定颗粒、三级颗粒膜、特定颗粒膜，主要涉及的分子功能为磷脂结合、酰胺结合、肽结合、整合素结合、β-淀粉样蛋白结合、细胞因子结合、货物受体活性、磷脂酰肌醇-3，4，5-三磷酸结合、肌肉结构成分、Toll样受体结合。KEGG通路富集分析结果显示，DEG主要涉及黑热病、破骨细胞分化、吞噬体、肺结核、血脂与动脉粥样硬化、冠状病毒、造血细胞系、细胞黏附分子、EB病毒感染、Rap1信号通路、金黄色葡萄球菌感染、肌动蛋白细胞骨架的调节、类风湿关节炎、白细胞跨内皮迁移、B细胞受体信号通路、溶酶体、自然杀伤细胞介导的细胞毒性、Fc epsilon RI信号通路、百日咳、炎症性肠病、抗原处理和提呈、补体和凝血级联反应、疟疾、军团菌病、病毒性心肌炎、过氧化物酶体增殖物激活受体（peroxisome proliferatorsactivated receptors，PPAR）信号通路、胆固醇代谢、哮喘、移植物抗宿主疾病、1型糖尿病等信号通路。

2.3 WGCNA结果

网络拓扑分析结果显示，软阈值为2，见图2。共识别到5个基因模块，其中绿松石色基因模块与AS相关性最强（r=-0.96，P＜0.001），见图3，该基因模块共包含15 087个基因。将绿松石色基因模块中的基因、FRG、DEG取交集，共得到17个交集基因，分别为ZEB1、ABCC5、CAPG、LGMN、CTSB、YAP1、PTPN6、BID、PLIN4、FTL、CDKN2A、PLIN2、HMOX1、IL1B、HAMP、ALOX5、DPP4，见图4。

图2 网络拓扑分析结果Figure 2 Network topology analysis results

图3 5个基因模块与AS的相关性Figure 3 Correlation between 5 gene modules and AS

图4 交集基因的韦恩图Figure 4 Venn diagram of intersecting genes

2.4 PPI分析和核心基因的确定

PPI分析结果显示，构建出1个含有17个节点、60条边的PPI网络图，见图5；核心基因分别为IL1B、CTSB、HMOX1、CDKN2A、ALOX5。

图5 交集基因的PPI网络图Figure 5 PPI network diagram of intersecting genes

2.5 核心基因验证

在GSE100927中，AS组IL1B、CTSB、HMOX1、CDKN2A、ALOX5表达水平高于对照组，差异有统计学意义（P＜0.05）；在GSE20129中，AS组IL1B、HMOX1表达水平高于对照组，差异有统计学意义（P＜0.05）；在GSE20129中，对照组与AS组CDKN2A、ALOX5表达水平比较，差异无统计学意义（P＞0.05）；在GSE226790中，对照组与AS组CTSB表达水平比较，差异无统计学意义（P＞0.05），见表1。

表1 不同数据集中对照组与AS组核心基因表达水平比较（±s）Table 1 Comparison of core genes expression levels between control group and AS group in different datasets

表1 不同数据集中对照组与AS组核心基因表达水平比较（±s）Table 1 Comparison of core genes expression levels between control group and AS group in different datasets

注：AS=动脉粥样硬化。

组别GSE100927GSE20129GSE226790样本量IL1BCTSBHMOX1CDKN2AALOX5样本量IL1BHMOX1CDKN2AALOX5样本量CTSB对照组126.10±0.98 8.83±0.67 8.43±1.00 6.17±0.51 7.77±1.02799.16±0.18 7.85±0.34 6.70±0.10 6.59±0.0935.87±1.03 AS组298.16±0.37 10.52±0.28 10.51±0.71 7.56±0.24 9.44±0.20569.50±0.29 8.30±0.19 6.73±0.10 6.64±0.0936.07±0.90 t值9.8711.467.5411.818.562.793.200.581.110.25 P值＜0.001＜0.001＜0.001＜0.001＜0.0010.0150.0060.5740.2870.817

2.6 ssGSEA结果

ssGSEA结果显示，IL1B、CTSB、HMOX1、CDKN2A、ALOX5主要涉及铁死亡、脂肪消化吸收、核酸结合寡聚结构域（nucleotide-binding oligomerization domain，NOD）样受体信号通路、组氨酸代谢、卟啉代谢、Th1和Th2细胞分化、TGF-β信号通路、FC epsilon RI信号通路等机制。

2.7 免疫浸润分析结果

对照组和AS组初始CD4T淋巴细胞、静息树突状细胞、活化自然杀伤细胞、嗜酸性粒细胞、初始B细胞表达水平均为0。AS组静息肥大细胞、静息自然杀伤细胞表达水平为0。AS组静息CD4记忆T淋巴细胞、浆细胞表达水平低于对照组，活化肥大细胞、单核细胞、滤泡辅助性T淋巴细胞、记忆B细胞表达水平高于对照组，差异有统计学意义（P＜0.05），见表2。

表2 对照组和AS组22种免疫浸润细胞表达水平比较（±s）Table 2 Comparison of expression levels of 22 types of immune infiltrating cells between the control group and the AS group

表2 对照组和AS组22种免疫浸润细胞表达水平比较（±s）Table 2 Comparison of expression levels of 22 types of immune infiltrating cells between the control group and the AS group

注：a表示U值；-表示无此项数据。

细胞对照组（n=12）AS组（n=29）t（U）值P值静息CD4记忆T淋巴细胞（images/BZ_14_2062_2150_2079_2187.png±s）0.148±0.0360.074±0.0305.197＜0.001活化CD4记忆T淋巴细胞（images/BZ_14_2062_2150_2079_2187.png±s）0.073±0.0530.042±0.0301.9690.058静息肥大细胞〔M（P25，P75）〕0.010（0，0.028）0——活化肥大细胞〔M（P25，P75）〕0.002（0，0.042）0.118（0.077，0.137）25.00a0.005活化树突状细胞〔M（P25，P75）〕0（0，0.078）0.011（0，0.030）82.00a0.828静息自然杀伤细胞〔M（P25，P75）〕0.008（0.002，0.036）0——中性粒细胞（images/BZ_14_2062_2150_2079_2187.png±s）0.083±0.0500.067±0.0310.9920.328单核细胞〔M（P25，P75）〕0.007（0，0.022）0.091（0.028，0.102）23.50a0.003浆细胞（images/BZ_14_2062_2150_2079_2187.png±s）0.011±0.0020.004±0.0033.7200.001 CD8 T淋巴细胞〔M（P25，P75）〕0.003（0，0.041）0.009（0，0.023）77.00a0.667调节性T淋巴细胞〔M（P25，P75）〕0（0，0.003）0（0，0.009）70.50a0.431滤泡辅助性T淋巴细胞〔M（P25，P75）〕0.009（0，0.030）0.031（0.020，0.045）37.00a0.027 γ-δT细胞（images/BZ_14_2062_2150_2079_2187.png±s）0.076±0.0380.074±0.0250.1420.888记忆B细胞（images/BZ_14_2062_2150_2079_2187.png±s）0.020±0.0080.033±0.0132.5450.016 M0型巨噬细胞（images/BZ_14_2062_2150_2079_2187.png±s）0.273±0.0420.302±0.0750.9160.366 M1型巨噬细胞〔M（P25，P75）〕0.024（0，0.037）0.007（0，0.016）52.00a0.126 M2型巨噬细胞（images/BZ_14_2062_2150_2079_2187.png±s）0.173±0.0880.134±0.0341.8490.073

相关性分析结果显示，IL1B表达水平与活化肥大细胞、中性粒细胞表达水平呈正相关，与滤泡辅助性T淋巴细胞、记忆B细胞、M0型巨噬细胞表达水平呈负相关（P＜0.05）；CTSB、HMOX1表达水平与22种免疫浸润细胞表达水平均无直线相关关系（P＞0.05）；CDKN2A表达水平与中性粒细胞表达水平呈负相关（P＜0.05）；ALOX5表达水平与活化CD4记忆T淋巴细胞、活化树突状细胞、M2型巨噬细胞表达水平呈负相关，与活化肥大细胞、调节性T淋巴细胞、M1型巨噬细胞表达水平呈正相关（P＜0.05），见表3。

表3 核心基因表达水平与22种免疫浸润细胞表达水平的相关性（r值）Table 3 Correlation between the expression level of core genes and the expression levels of 22 types of immune infiltrating cells

3 讨论

自从DIXON等［13］在2012年首次提出铁死亡概念以来，铁死亡与各种疾病之间的关系被广泛研究［14-15］。研究表明，铁死亡在心血管疾病中起关键作用，包括AS［16］、心肌梗死［17］和心力衰竭［18］等。但哪些FRG与AS有关，目前尚不清楚。本研究旨在基于WGCNA筛选AS中铁死亡核心基因，并分析其与免疫浸润细胞的关系。

本研究通过cytoHubba插件的Degree算法筛选出了5个AS中铁死亡核心基因，分别为IL1B、CTSB、HMOX1、CDKN2A、ALOX5。进一步验证核心基因，结果显示，在GSE100927中，AS组IL1B、CTSB、HMOX1、CDKN2A、ALOX5表达水平高于对照组；在GSE20129中，AS组IL1B、HMOX1表达水平高于对照组；在GSE20129中，对照组与AS组CDKN2A、ALOX5表达水平比较，差异无统计学意义；在GSE226790中，对照组与AS组CTSB表达水平比较，差异无统计学意义；分析原因，可能与所用的数据集样本量较小、样本来源不同以及样本间存在差异等多种因素有关。本研究ssGSEA结果显示，IL1B、CTSB、HMOX1、CDKN2A、ALOX5主要涉及铁死亡、脂肪消化吸收等作用，说明核心基因通过多因素、多机制共同参与AS的发生发展过程。研究显示，白介素1β（interleukin-1β，IL-1β）由单核细胞、星形胶质细胞和脑内皮细胞释放［19］，其对内皮细胞、单核细胞以及血管平滑肌细胞均有促进作用，可促进AS的发生。IL-1β还可通过刺激冠状动脉壁内皮细胞而产生细胞间黏附分子和血管黏附分子，从而促进单核细胞趋化蛋白的合成，进而导致单核细胞的聚集，进一步推动AS的进展［20］。而铁死亡可通过抑制GPX4和SIRT1的表达来增加IL1B的表达，进而激活炎症反应，产生大量活性氧，最终导致AS的发生［21］。CTSB是一种广泛表达于溶酶体上的半胱氨酸肽链内切酶，在自噬、抗原提呈、细胞应激、新陈代谢和溶酶体依赖性细胞死亡等方面发挥着不可或缺的作用［22］，其可通过加重炎症反应和降解含弹性蛋白的动脉细胞外基质而参与动脉的变性、钙化和硬化［23］。研究显示，在冠状动脉病变组织中，与自噬相关的CTSB的表达明显上调［24］；CTSB在人颈动脉斑块中也呈高表达，并与斑块严重程度相关［25］。此外，血糖升高和肥胖是公认的AS的危险因素，而CTSB在肥胖、糖尿病、非酒精性脂肪性肝病和癌症等肠道血管生成和胆固醇吸收中起关键作用［26］。HMOX1是血红素降解过程的限速酶，可控制一氧化碳、Fe2+和胆绿素的生成，其与AS的发生发展关系密切，可抑制各种炎症因子导致的内皮细胞凋亡，从而抑制AS形成，并可促进晚期斑块的稳定［27］。CDKN2A可参与细胞增殖、凋亡及调节细胞周期［28］。p16INK4a、p14ARF和p15INK4b为CDKN2A/B基因位点的蛋白质产物，在AS病变中，CD68阳性巨噬细胞和平滑肌细胞可分泌这些蛋白质产物，而其均可参与AS的形成［29］。ALOX5是一种含铁的非血红素双加氧酶，可激活炎症反应并触发各种细胞死亡模式［30-31］。已有研究证实，ALOX5是引发癌症［32］或神经元损伤［33］中铁死亡的关键酶。但尚未见有直接证据表明ALOX5可介导AS斑块中的铁死亡，这仍需要进一步研究证实。

免疫浸润细胞在AS中持续活化和分化，这对细胞间抗原免疫提呈和细胞信号通路的激活具有重要作用［34］，因而了解免疫微环境特征有利于认识AS斑块生物学过程。本研究免疫浸润分析结果显示，AS组静息CD4记忆T淋巴细胞、浆细胞表达水平低于对照组，活化肥大细胞、单核细胞、滤泡辅助性T淋巴细胞、记忆B细胞表达水平高于对照组，提示AS患者免疫微环境发生了改变。研究表明，IL1B可导致AS患者体内中性粒细胞明显增加［35］，并促进滤泡辅助性T淋巴细胞的分化［36］；CDKN2A可阻止中性粒细胞胞外诱捕网（neutrophil extracellular traps，NETs）的产生，从而减缓AS的发生发展［37-38］；ALOX5可通过促进M1型巨噬细胞的产生而参与AS斑块的形成［39-40］，还可通过抑制M2型巨噬细胞的生成而发挥抗炎作用，从而减缓AS的病程进展［40-41］。本研究相关性分析结果显示，IL1B表达水平与活化肥大细胞、中性粒细胞表达水平呈正相关，与滤泡辅助性T淋巴细胞、记忆B细胞、M0型巨噬细胞表达水平呈负相关；CTSB、HMOX1表达水平与22种免疫浸润细胞表达水平均无直线相关关系；CDKN2A表达水平与中性粒细胞表达水平呈负相关；ALOX5表达水平与活化CD4记忆T淋巴细胞、活化树突状细胞、M2型巨噬细胞表达水平呈负相关，与活化肥大细胞、调节性T淋巴细胞、M1型巨噬细胞表达水平呈正相关；与既往研究结果［42-43］存在差异，可能与样本量较小有关，后续仍需进行更全面的数据分析及验证。

4 结论

综上所述，本研究基于WGCNA共筛选出5个AS中铁死亡相关核心基因，分别为IL1B、CTSB、HMOX1、CDKN2A、ALOX5，其中IL1B、CDKN2A、ALOX5表达水平与部分免疫浸润细胞表达水平有直线相关关系，而CTSB、HMOX1表达水平与22种免疫浸润细胞表达水平均无直线相关关系。但本研究未能纳入AS患者的一些临床特点进行综合分析，如年龄、性别和饮食习惯等，且未通过蛋白表达的检测来验证核心基因的翻译水平。

作者贡献：方柔柔、赵晶提出主要研究目标，负责研究的构思与设计；方柔柔进行研究的实施，撰写及修订论文；方柔柔、杨启帆、韩若冰、孙娜进行数据的收集与整理，统计学处理，图、表的绘制与展示；方柔柔、邬东东、李娟、徐守竹进行结果的分析与解释；赵晶负责文章的质量控制与审查，对文章整体负责、监督管理。

本文无利益冲突。