基于加權基因共表達網絡分析的動脈粥樣硬化中鐵死亡核心基因及其與免疫浸潤細胞的關系研究

方柔柔，楊啟帆，韓若冰，鄔東東，孫娜，李娟，徐守竹，趙晶

作者單位：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患者的一些臨床特點進行綜合分析，如年齡、性別和飲食習慣等，且未通過蛋白表達的檢測來驗證核心基因的翻譯水平。

作者貢獻：方柔柔、趙晶提出主要研究目標，負責研究的構思與設計；方柔柔進行研究的實施，撰寫及修訂論文；方柔柔、楊啟帆、韓若冰、孫娜進行數據的收集與整理，統計學處理，圖、表的繪制與展示；方柔柔、鄔東東、李娟、徐守竹進行結果的分析與解釋；趙晶負責文章的質量控制與審查，對文章整體負責、監督管理。

本文無利益沖突。