动脉粥样硬化斑块破裂的关键免疫相关基因和潜在治疗药物的识别

尚莎莎，陈玉善，李晓辉，王建茹

（河南中医药大学 第一附属医院 心脏中心，河南 郑州 450000）

动脉粥样硬化（atherosclerosis，AS）是一种涉及免疫反应的慢性炎症性疾病，可引起许多临床并发症，如冠状动脉疾病、外周动脉疾病和缺血性脑卒中［1-3］。近几十年来，能够有效预防或抑制AS病理过程的临床治疗方法有了长足的进展。然而，由动脉粥样硬化斑块破裂引起的急性心肌梗死和缺血性卒中仍是世界上相当大的死亡原因［4-6］。尽管，诸多研究者对AS的病理机制进行了大量的探索研究，但对其仍知之甚少，尤其是斑块破裂的机制。免疫系统在AS的发生发展过程中起着重要作用，是研究防治AS的重要切入点［3］。然而，参与动脉粥样硬化斑块破裂（atherosclerotic plaque rupture，APR）的特异性免疫机制尚不完全明确，斑块破裂的生物学标记物也缺乏。在此背景下，研究APR病理进展的免疫机制，寻找合适的靶点防治斑块破裂引起的临床并发症具有重要意义。

现阶段，已有诸多学者对人类AS不同阶段和病变部位的表达谱进行了研究，包括颈动脉AS斑块（GSE24495、GSE118481和GSE43292）、外周血管AS斑块（GSE23304、GSE24702）、冠状动脉AS 斑 块（GSE11138）、颈动脉钙化斑块（GSE125771，GSE159677）、颈动脉早期和晚期AS斑块（GSE28829）、颈动脉晚期AS斑块（GSE97210）、颈动脉破裂和稳定 AS斑块（GSE120521，GSE41571和E-MTAB-2055）等。同时，许多研究也通过将基因芯片和高通量测序数据与生物信息学分析相结合的方法，来探索AS可能的病理进展机制，并取得了一些进展［7-16］。然而，有关APR的潜在免疫机制的研究相对较少。因此，积极探索稳定动脉粥样硬化斑块进展为破裂斑块的免疫景观，挖掘能够早期诊断ARP的免疫相关诊断标志物非常重要。

本研究从基因表达综合数据库（gene expression omnibus，GEO）和ArrayExpress数据库中下载与稳定动脉粥样硬化斑块（stable atherosclerotic plaque，SAP）和破裂动脉粥样硬化斑块（ruptured atherosclerotic plaque，RAP）相关的数据集（GSE120521、GSE41571和E-MTAB-2055）。利用SVA包将3个数据集合并为一个数据集，并利用Limma包筛选合并数据集中SAP和RAP样本间的差异表达基因（differentially expressed genes，DEGs），再与免疫相关基因（immune-related genes，IRGs）取交集，获取差异表达的免疫相关基因（differentially expressed immune-related gene，DEIRGs）。然 后，构 建DEIRGs的蛋白-蛋白相互作用（protein-protein interaction，PPI）网络，再利用分子复合物检测（molecular complex detection，MCODE）插件筛选hub基因。利用连接图谱（connectivity Map，cMap）预测了靶向DEIRGs，治疗AS斑块破裂的候选化合物，并将其与hub基因进行了分子对接，来评价他们之间的作用关系。同时，评价了hub基因的诊断效能，并开展了hub基因的单基因基因集富集分析（gene set enrichment analysis，GSEA）。本研究的思路流程见图1。本研究将有助于了解SAP和RAP病理过程中潜在的生物标志物及免疫分子机制，指导更精确的靶向免疫治疗，从而进一步为其防治提供新的思路和科学依据。

图1 研究思路流程图Figure 1 Flowchart of the study design

1 材料与方法

1.1 数据的来源和预处理

为了获得与AS进展相关的基因，从GEO数据 库（http：／／www.ncbi.nlm.nih.gov／geo）和ArrayExpress数据库（https：／／www.ebi.ac.uk／arrayexpress／）中下载由SAP和破裂动脉粥样硬化斑块组成的数据集。GSE120521，GSE41571和E-MTAB-2055 3个数据集被采用，AS斑块来自人颈动脉。3个数据集的详细信息详见表1。分别下载3个数据集的矩阵文件和平台注释文件。将GSE120521数据集的FPKM数据转换为TPM数据。根据平台注释文件，通过Perl脚本将GSE41571和E-MTAB-2055数据集中的探针名称转换为相应的基因名。如果多个探针对应一个基因，则这些探针的平均值即是该基因的最终表达量。将3个数据集合并为一个数据集（即合并数据集），并使用SVA包消除批次效应。

表1 3个选定的数据集的相关信息Table 1 The relevant information of three selected datasets

1.2 DEGs基因的鉴定及其功能富集分析

利用limma包，以｜log2（fold change，FC）｜＞1和P＜0.05为阈值标准对合并数据集的表达数据进行差异分析，筛选DEGs［17］。然后，再利用clusterProfiler包对DEGs进行基因本体论-生物过程（gene ontology，GO）和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析。P＜0.05被认为是显著的富集。

1.3 DEIRGs的鉴定

分别从import数据库（https：／／www.immport.org／shared／genelists）和MsigDB数据库（http：／／www.gsea-msigdb.org／gsea／index.jsp）提取IRGs。将筛选出的DEGs和获取的IRGs取交集，获取DEIRGs。

1.4 DEIRGs的蛋白互作网络的构建、模块分析及hub基因的筛选

将获得的DEIRGs上传至STRING数据库（https：／／string-db.org／），选择“Multiple proteins”模式，互作得分设置为最高置信度≥0.400，进行PPI分析，并将结果导入Cytoscape3.7.2软件构建PPI网络。利用MCODE插件并以其默认参数（degree cut-off＝2，node score cut-off＝0.2，Max depth＝100，k-score＝2）对PPI网络进行Module分析，筛选最重要的子模块，最重要子模块中的基因被认为是hub基因。利用Pearson相关分析对hub基因间的相关性进行分析。此外，利用pROC包绘制ROC曲线，并计算曲线下面积（area under the curves，AUC），来评价各hub基因的诊断效能。

1.5 cMap分析和分子对接

依据GPL96平台（Affymetrix Human Genome U133A Array），通过Perl脚本将DEIRGs转换成对应的特征集ID，然后将这些ID上传到cMap数据 库（https：／／portals.broadinstitute.org／cmap／）中。连接性分数（范围是-1到+1）来评价化合物和DEIRGs之间的相似性。连接性分数负值越高的化合物代表其对AS进展的潜在治疗作用越大。筛选阈值设定为P＜0.01，n≥3，均值＜-0.3。从PubChem（https：／／pubchem.ncbi.nlm.nih.gov／）数据库中下载化合物的2D结构和分子式；从PDB（http：／／www.rcsb.org／）数据库中下载hub基因的3D结构。最后，利用Autodock Vina软件进行分子对接，获取结合能值，对接结果采用Discovery Studio Visualizer 4.5软件进行可视化。结合能数值为负则提示活性成分与hub基因可自由结合，数值越小说明两者的作用关系越可靠。

1.6 单个hub基因的GSEA

以单个hub基因在合并数据集所有样本表达量的中位数为依据，将数据集中的基因分为高表达组和低表达组。然后，利用GSEA4.1.0软件，将排列的数目设置为1 000，选择c2.cp.kegg.v7.4.symbols.gmt作为参考基因集进行GSEA。最后，按照｜NES｜＞1、P＜0.05和FDR＜0.25为阈值筛选阳性基因集。

2 结果

2.1 DEGs的鉴定及其功能富集分析

由图2可知SVA包合并后数据间的批次效应已不是很明显。共鉴定出93个下调和98个上调基因。GO分析结果显示，依据P＜0.05确定了191个GO条目，其中DEGs在嗜中性粒细胞脱粒、嗜中性粒细胞活化参与免疫应答、嗜中性粒细胞介导免疫、抗原加工和呈递等过程中显著富集（图3A）。KEGG通路富集分析结果显示，依据P＜0.05共映射出了35条信号通路，主要涉及溶酶体、吞噬体、胆固醇代谢、抗原加工和呈递、补体和凝血级联、细胞黏附分子、自然杀伤细胞介导的细胞毒性等（图3B）。

图2 3个数据集合并前后的结果Figure 2 The results of three datasets before and after merging

图3 DEGs的功能富集分析的结果Figure 3 The results of functional enrichment analyses of DEGs

2.2 DEIRGs的鉴定

从MSigDB数据库中的获得Immune_response（M19817）和Immune_system_process（M13664）两个免疫相关基因集合，分别包含235和332个IRGs；从ImmPort数据库中获得1793个IRGs；合并并剔除重复后共获得1959个IRGs。将191个DEGs和1959个IRGs取交集后，共获得29个DEIRGs，其中上调基因26个，下调基因3个（图4）。

图4 DEIRGs的韦恩图Figure 4 Venn diagram of DEIRGs.

2.3 cMap分析

将DEIRGs上传到cMap数据库，结果共获得1309个化合物，其中有9个候选化合物满足筛选的条件（表2）。9个候选化合物的化学结构见图5。

图5 候选化合物的2D结构Figure 5 2Dstructure of the candidate compounds

表2 通过cMap数据库鉴定出的潜在化合物Table 2 The potential compounds identified by the cMap database.

2.4 PPI网络构建、模块分析及hub基因的鉴定

29个DEIRGs所构建PPI网络包括72条边和27个节点（图6A）。利用MCODE插件对PPI网络进行聚类并构建功能模块，结果共筛选出1个重要的子模块，该模块包括8个节点和25条边（图6B）。这8个节点即是hub基因，包括单核细胞分化抗原CD14（monocyte differentiation antigen CD14，CD14）、集落刺激因子1受体（colonystimulating factor 1 receptor，CSF1R）、高亲和力免疫球蛋白受体γ亚型（high-affinity immunoglobulin epsilon receptor subunit gamma，FCER1G）、S100钙结合蛋白A9（s100 calcium-binding protein A9，S100A9）、低亲和力免疫球蛋白γFc区受体iii-A（low-affinity immunoglobulin gamma Fc region receptorⅢ-A，FCGR3A）、造血细胞信号转换器（hematopoietic cell signal transducer，HCST）、TYRO蛋白酪氨酸激酶结合蛋白（TYRO protein tyrosine kinase binding protein，TYROBP）、整合素亚基β-2（integrin beta 2，ITGB2）。所有hub基因均呈正相关，其中ITGB2与TYROBP的相关性最强（图6C）。AUCs结果表明，所有hub基因在区分SAP和RAP方面均表现出良好的诊断性能，可作为AS斑块破裂的潜在诊断生物标志物（图6D）。

图6 hub基因的鉴定Figure 6 Identification of hub genes

2.5 分子对接的结果

为了验证cMap的分析结果，将hub基因与其对应的潜在化合物进行分子对接。图7A所示，共有72对“潜在化合物-作用靶点”关系对，且结合自由能均小于-5 kcal／moL，这提示筛选出来的hub基因与其对应的潜在化合物具有良好的亲和力，表明cMap分析结果具有一定的可信度。图7B-E展示了结合能小于-8 kcal／moL的对接结果。

图7 hub基因与其对应潜在化合物分子对接的3D和2D示意图Figure 7 Schematic（3D and 2D）representation that molecular docking of the hub genes with its corresponding components

2.6 单基因GSEA结果

为了进一步阐释hub基因在AS斑块破裂进展过程中可能的作用机制，本研究在合并数据集中分别以单个hub基因开展了GSEA。研究结果显示，8个hub基因在其低表达组中均未富集到阳性基因集；在CD14、CSF1R、FCER1G、FCGR3A、HCST、ITGB2、S100A9和TYROBP的高表达组中分别富集到了42、42、34、43、44、39、43和38个阳性基因集，取交集后共有16个共同基因集，主要包括自然杀伤细胞介导的细胞毒性、抗原加工和呈递、γ介导的吞噬作用、B细胞受体信号通路、趋化因子信号通路、NOD样受体信号通等（图8）。

图8 单个hub基因GSEA的upset图Figure 8 The upset plot of single-gene GSEA of hub genes

3 讨论

AS是冠状动脉疾病、缺血性中风等许多动脉粥样硬化性疾病的共同病理生理机制的基础，这些疾病是世界上相当多的死亡和发病的主要原因。AS的发生发展是一个多因素参与、多步骤的复杂病理过程。既往大量研究表明，先天免疫和适应性免疫机制均可加速或抑制AS，其已成为近年来研究AS病理过程的热点［1、3、18］。因此，在这一研究领域应深入探索，以便发现相应的新的防治措施。

为了避免单数据集分析存在个体偏倚的情况，本研究使用SVA包将3个符合条件的数据集合并为了1个数据集，并在合并数据集中鉴定出191个DEGs。GO和KEGG富集分析结果表明，DEGs主要参与了免疫应答、抗原加工递呈、细胞黏附和自然杀伤细胞介导的细胞毒性等通路和功能。研究显示，中心粒细胞的活化脱颗粒及其介导的免疫反应、抗原加工和呈递等炎症免疫机制参与了AS的病理过程［19-20］，这些结果与本次研究的结果相似。以上结果提示，DEGs调控AS恶化的机制与免疫相关的生物学过程和通路有关。

众所周知，与APR相关的hub基因的识别对严重AS疾病的诊断和防治非常重要。本研究共发现了29个DEIRGs，从中筛选出了8个hub基因。AUCs结果表明，8个hub基因在区分SAP和RAP方面均表现出了良好的诊断性能，是APR的候选诊断生物标志物。经查阅文献，8个hub基因中只有CD14和S100A9被报道参与了APR。Hermansson等［21］报道CD14在血栓性颈动脉斑块中的表达显著增加。已有研究表明，S100A9（又称MRP-14）被鉴定为易破裂动脉粥样硬化斑块的局部标记物；S100A8／A9复合体是急性心肌梗死后局部炎症反应的生物标记物，与冠状动脉和颈动脉粥样硬化程度具有一定的相关［22-24］。此外，HCST在AS病理过程中的作用机制尚未见相关报道。虽然，已有文献报道了其他hub基因（CSF1R［25］、FCER1G［26］，FCGR3A［27］、ITGB2［28-29］和TYROBP［29］）介导了AS的病理过程。但是，它们在斑块破裂中的致病作用机制尚未被报道。

为了更好地研究hub基因在AS恶化中的潜在生物学功能，本研究对单个hub基因分别开展了GSEA。GSEA结果表明，hub基因主要参与了免疫相关通路，包括自然杀伤（NK）细胞介导的细胞毒性、抗原加工和呈递、细胞黏附分子、B细胞受体（BCR）信号通路、核苷酸结合寡聚域（NOD）样受体信号通路等。回顾文献发现上述提及的一些通路被证实参与了AS的病理过程。NOD1和NOD2是NOD样受体家族的创始成员，被认为是AS因子，尤其是NOD1在APR方面发挥了突出作用［30-31］。研究表明，NK细胞介导的细胞毒性功能异常与AS心血管疾病密切相关［32］。NK细胞可通过细胞毒性依赖机制促进AS和AS坏死核的扩张［33］。B细胞通过细胞间接触、抗原呈递和细胞因子的产生来调节细胞免疫反应，进而参与促AS病理过程中的免疫反应［34］。此外，在AS病理过程中，BCR依赖的信号通过识别氧化特异性抗原表位的相关抗原而被激活［34］。总之，通过文献复习，不难发现这些hub基因富集出的通路与既往文献报道的结果存在一些吻合，在一定程度上佐证了本研究结果的科学性和准确性。

随之，本研究利用cMap数据库筛选了潜在延缓逆转稳定斑块和破裂斑块之间DEIRGs异常表达的小分子，并获得了9个候选化合物，包括阿普萘洛尔（alprenolol）、氟伏沙明（fluvoxamine）、噻嘧啶（pyrantel）、石蒜碱（lycorine）、硫马唑（sulmazole）、苯海索（trihexyphenidyl）、加贝酯（gabexate）、链霉素（streptomycin）、吗啉胍（moroxydine）。氧化低密度脂蛋白与不稳定斑块坏死核的发生和进展具有很强的相关性，阿普萘洛尔虽然是一种β受体阻滞剂，也可抑制低密度脂蛋白的氧化修饰，达到延缓或抑制AS恶化的作用［35］。氟伏沙明作为一种选择性5-羟色胺再吸收抑制剂型的抗抑郁药，临床上已被广泛用于各种抑郁症的治疗，但其抗AS的确切分子机制尚不清楚。细菌、寄生虫、病毒等病原体感染参与了AS的发生发展，为探索AS的病理机制及干预措施提供了研究思路［36-37］。基于上述认识，吗啉胍（广谱抗病毒药）、噻嘧啶（广谱驱虫药）、链霉素（氨基糖苷类抗生素）虽目前尚未有文献报道其可防治AS，但表现出了潜在的干预作用。石蒜碱是药用石蒜科植物中的重要异喹啉类生物碱，具有抗肿瘤、抗病毒、抗炎、调控脂质合成、抗寄生虫等多种生物学活性［38］。研究显示，石蒜碱对心血管疾病也表现出了一定的防治潜力，其不仅可以在体内外实验中有效改善心肌纤维化，还可通过抑制炎症、氧化应激和细胞凋亡等来改善异丙肾上腺素诱导的心功能不全［39-40］。加贝酯是一种非肽类的蛋白酶抑制剂，临床常用于胰腺炎、弥漫性血管内凝血等，其在心血管方面的应用主要为促进心脏缺血后功能的恢复［41］，左心室辅助装置安装后的抗凝治疗［42］。硫马唑作为苯咪唑类衍生物磷酸二酯酶抑制剂，具有明显的正性肌力和血管扩张作用，并能增强心肌对Ca2+的敏感性，可用于充血性心衰等疾病。研究报道，硫马唑能改善缺血再灌注后顿抑心肌的收缩力及能量代谢，减少心律失常及心肌细胞损伤［43］。苯海索作为一种中枢抗胆碱类药物，临床常用于帕金森病、帕金森综合征及药物引起的锥体外系疾病。早期国内有研究报道，苯海索可对抗乌头碱引起的室性心律失常，还可抑制心肌Na+、Ca2+跨膜转运作用［44-45］。早期国外学者曾报道苯海索成功治疗严重变异性心绞痛的案例［46］。迄今为止，除了阿普萘洛尔外，其他候选化合物尚未见被报道用于AS的治疗，但本研究的结果可为这些“老药物”拓展“新用途”提供了一些应用思路。此外，本研究分子对接结果显示，候选化合物与相应hub基因之间具有良好的结合能力，这从侧面佐证了cMap分析结果的科学性和可信性。综上所述，这些候选化合物表现出了潜在的抑制APR的作用，未来无疑值得深入研究。

本研究尚有一些的局限性，主要体现在以下几个方面：①本研究纳入分析的样本量有限。尽管本研究使用SVA包合并了3个数据集来弥补这个限制，以保证结果的客观性；但因为开放的公共数据库中满足研究条件的数据集较少，因此研究结果仍需进一步在大样本量的数据中进行验证。②cMap数据库虽然已被广泛用于用于揭示小分子化合物、基因以及疾病状态的功能联系，但仍跟实际情况存在一些偏倚。③本研究没有对结果进行实验验证，因此仍需要后续的基础和临床研究来进行佐证。

总之，本研究利用生物信息学分析和分子对接技术，鉴定出了8个hub基因在APR的进展过程中起关键作用，其可能是APR的潜在诊断生物标志物；发现了9个候选化合物表现出了靶向免疫机制治疗APR的潜在作用。本研究的核心价值和意义主要体现在以下几个方面：①APR的特异性免疫机制尚不完全明确，本研究有助于更好地从免疫学角度来探索APR的病理机制；②目前，用于诊断APR的生物学标记物尚缺乏，早期或及时的诊断对于APR的防治尤为重要，本研究可为APR潜在临床诊断标志物的开发提供了一定的借鉴；③靶向免疫机制来治疗AS一直是心血管领域研究的重点和热点，本研究所挖掘出的DEIRGs，尤其是hub基因为靶向免疫机制防治APR提供了方向；④本研究所挖掘出的候选化合物，为潜在抑制APR药物的开发提供了一定的科学依据，同时也为“老药新用”的应用提供了新的方向和思路。

作者贡献声明

尚莎莎：数据的整理、文章撰写；陈玉善：提出研究思路，论文审阅及修订；李晓辉：负责数据集的下载、cMap分析和分子对接的工作；王建茹：提出研究思路，完成数据的校正、差异基因及关键基因的筛选、功能富集分析等生物信息学分析工作。

利益冲突声明

本研究未受到企业、公司等第三方资助，不存在潜在利益冲突。