基于生物信息学分析类风湿关节炎与急性心肌梗死共同差异基因及相关机制

Analysis of Common Differential Genes and Related Mechanisms of Rheumatoid Arthritis and Acute Myocardial Infarction Based on Bioinformatics

摘要: 目的：基于生物信息学分析类风湿关节炎(RA)与急性心肌梗死(AMI)共同差异基因、共同通路、免疫机制及潜在药物。方法：通过GEO数据库下载RA数据集(GSE77298)及AMI数据集(GSE66360)，应用R语言分别筛选RA及AMI数据集的差异表达基因(DEGs)。对DEGs进行GO和KEGG富集分析。利用韦恩图对RA和AMI差异表达基因取交集筛选共同DEGs，利用CytoScape软件构建共同DEGs的PPI网络识别参与疾病发生发展的关键基因(HubGene)。HubGene在外部数据集验证。对HubGene进行KEGG、GSEA富集分析并导入TRRUST数据库预测出转录调控因子，对HubGene行ssGSEA评估与免疫细胞关系。将关键基因导入DGIdb数据库寻找潜在治疗药物。结果：RA数据集共鉴定出250个DEGs，富集分析主要集中在免疫、趋化因子活性、细胞分泌等。AMI数据集共鉴定出438个DEGs，富集分析主要集中在免疫、炎症、感染等相关生物学改变。通过CytoScape筛选出关键基因IL1RN，对IL1RN进行单基因KEGG、GSEA富集分析，均富集到NF-κB信号通路，在TRRUST数据库预测出转录调控因子NFκB1，ssGSEA结果提示IL1RN与多种免疫细胞相关，DGIdb数据库得到甲氨蝶呤、双醋瑞因、低分子肝素、氟哌啶醇四种药物。结论：IL1RN可能是AS发生AMI的预测指标，而IL1RN与NFκB1存在负相关，进而影响NF-κB信号通路表达。并且IL1RN的表达与免疫细胞呈相关性，甲氨蝶呤、双醋瑞因、低分子肝素、氟哌啶醇可能是AS合并AMI的潜在药物。

Abstract: Objective: To analyze the common differential genes, common pathways, immune mechanisms and potential drugs between rheumatoid arthritis and acute myocardial infarction based on bioinformatics. Methods: RA dataset (GSE77298) and AMI dataset (GSE66360) were downloaded from the GEO database, and R language was used to screen differentially expressed genes (DEGs) in RA and AMI datasets. GO and KEGG enrichment analyses were performed for DEGs. Venn diagram was used to screen common DEGs between RA and AMI differentially expressed genes, and CytoScape software was used to construct PPI network of common DEGs to identify the HubGene involved in the occurrence and development of the disease. HubGene was validated on an external dataset. KEGG and GSEA enrichment analysis were performed on HubGene and imported into TRRUST database to predict transcription regulators. ssGSEA was used to evaluate the relationship between HubGene and immune cells. Key genes were imported into DGIdb database to find potential therapeutic drugs. Results: A total of 250 DEGs were identified in the RA dataset, and enrichment analysis mainly focused on immunity, chemokine activity, and cell secretion. A total of 438 DEGs were identified in the AMI dataset, and enrichment analysis mainly focused on immune, inflammation, infection and other related biological changes. The key gene IL1RN was screened out by CytoScape, and the KEGG and GSEA enrichment analysis of IL1RN showed that all of them were enriched in the NF-κB signaling pathway. The transcription regulator NFκB1 was predicted in TRRUST database. Four drugs were obtained from the drug-drug-gene interaction database (DGIdb): methotrexate, diacerein, tinzaparin sodium and haloperidol decanoate. Conclusion: IL1RN may be a predictor of AMI in AS patients, and IL1RN is negatively correlated with NFκB1, thereby affecting the expression of NF-κB signaling pathway. The expression of IL1RN is correlated with immune cells. Methotrexate, diacerein, tinzaparin sodium and haloperidol decanoate may be potential drugs for AS complicated with AMI.

1. 引言

据报道，慢性全身性炎症性疾病(CSIDs)如类风湿性关节炎(rheumatoid arthritis, RA)与包括急性心肌梗死(acute myocardial infarction, AMI)在内的缺血性心血管事件的风险增加相关 [1] [2] 。慢性炎症及其相关机制包括髓内功能障碍、氧化应激、巨噬细胞积累和促炎细胞因子被认为是CSIDs [3] [4] [5] 心血管疾病发病机制的关键特征。一篇meta分析 [6] 纳入24项观察性研究、共111,758例患者，结论是RA患者死于粥样硬化性冠状动脉疾病(coronary artery disease, CAD)的风险比一般人群高59%。当患者发生包括急性心肌梗死在内的冠状动脉疾病时，有CSIDs的患者的预后比没有的患者更差 [7] 。随着微阵列技术的不断发展和进步，可通过研究差异表达基因(DEGs)来筛选疾病的候选生物标志物，为疾病的治疗方法开拓新思路 [8] 。在这项研究中，笔者通过生物信息学分析和微阵列技术，研究GEO数据库中RA与AMI相关的数据集，探索RA与AMI的中枢基因可能为疾病的诊断治疗提供新的生物信息学依据，获得的潜在生物标志物和信号通路可能有助于解释RA与AMI之间的联系。

2. 资料与方法

2.1. 数据收集

从基因表达数据库(Gene expression omnibus, GEO, https://www.ncbi.nlm.nih.gov/geo/) [9] 下载RA基因表达数据集GSE77298、GSE55235作为测试集和验证集，GSE77298包含16例RA滑膜组织和7例正常人滑膜组织，GSE55235包含10例RA滑膜组织和10例正常人滑膜组织；下载AMI基因表达数据集GSE66360、GSE48060作为测试集和验证集，GSE66360包含49例AMI患者外周血分离出的循环内皮细胞和50例正常人外周血分离出的循环内皮细胞，GSE48060包含31例AMI患者外周血样本和21例正常人外周血样本，GSE77298、GSE66360、GSE48060均来自于GPL570芯片平台，GSE55235来自于GPL96芯片平台。

2.2. DEGs筛选

利用R语言软件中的“limma”包对GSE77298和GSE66360数据集进行标准化处理，再用GPL570平台将探针转换成基因名称，以|logFC (fold change)| ≥ 1和调整后的P < 0.05为筛选标准，筛选DEGs。结果用R语言“ggplot”和“pheatmap”包绘制火山图和热图展示。用维恩图展示他们的共同DEGs。

2.3. DEGs功能富集分析

利用R软件“clusterProfiler”包 [10] 对DEGs进行在基因本体论(Gene Ontology, GO)和京都基因和基因组百科全书(Kyoto Encyclopedia of Genes and Genomes, KEGG)功能富集分析。GO分析分别从生物过程(Biological process, BP)、细胞成分(Cellular component, CC)和分子功能(Molecular function, MF)三个部分进行富集，KEGG分析对DEGs进行信号通路富集，分别对排行的前5位绘制柱状图。

2.4. 蛋白质与蛋白质相互作用网络(Protein-Protein Interaction Networks, PPI)构建

将筛选出来的DEGs导入STRING网站(https://www.string-db.org) [11] 以综合得分超过0.4作为互作的筛选条件进行分析获得PPI网络。利用Cytoscape [12] 软件中的Degree (连接度)算法对PPI网络进行关键基因(HubGene)鉴定。

2.5. HubGene外部数据集验证

筛选出的HubGene在GSE55235、GSE48060数据集验证。以调整后的P < 0.05为筛选条件。运用R语言软件中“ggpubr”包进行可视化处理。

2.6. 预测信号通路和转录调控因子

以HubGene表达情况进行KEGG富集分析及基因集富集分析(Gene Set Enrichment Analysis, GSEA)并导入转录因子调控网络数据库(Trrust) (https://www.grnpedia.org/trrust/)预测出调控核心基因的转录因子(transcription factor, TF) [13] 。

2.7. HubGene与免疫细胞浸润水平的相关性分析

以HubGene表达情况使用R语言“GSVA包”进行28种免疫细胞的单样本基因集富集分析(ssGSEA) [14] ，评估免疫浸润程度。

2.8. 预测潜在药物

DGIdb在线数据库(https://dgidb.org/)整合了现有30多个数据库中药物与基因互作关系，将得到的差异基因输入DGIdb数据库挖掘HubGene潜在相互作用药物。

3. 结果

3.1. RA和AMI差异基因获取

RA数据集共鉴定出250个DEGs，其中包含163个上调基因和87个下调基因。AMI数据集共鉴定出438个DEGs，其中包含344个上调基因和94个下调基因。使用维恩图鉴定出两种疾病的14个交集基因，分别是IL1RN、CCR5、TNFAIP6、GCSAM、CRTAM、ZNF382、MEST、TREM1、BCL2A1、RGS1、KCNJ15、PTX3、C15orf48、LILRB3 (图1(a)~(e))。

3.2. DEGs功能富集分析

对RA数据集DEGs进行GO和KEGG功能富集分析；GO富集分析中与BP相关的主要集中在基于免疫球蛋白超家族结构域构建的免疫受体体细胞重组的适应性免疫应答(adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains)、肌浆网钙离子转运(sarcoplasmic reticulum calcium ion transport)、出牙(tooth eruption)等。与CC相关的主要集中在肌原纤维(myofibril)、可收缩纤维(contractile fiber)、肌原纤维节(sarcomere)等。与MF相关的主要集中硫酸腺苷酸激酶活性(adenylylsulfate kinase activity)、硫酸腺苷基转移酶活性(sulfate adenylyltransferase activity)、硫酸腺苷酰转移酶(ATP)活性(sulfate adenylyltransferase (ATP) activity)等。KEGG主要富集的通路是胰岛素分泌(Insulin secretion)、胃酸分泌(Gastric acid secretion)、脂肪细胞中脂肪分解的调节(Regulation of lipolysis in adipocytes)等(图2(a)，图2(b))。

(a) (b) (c) (d) (e)

Figure 1. (a) Volcano map of DEGs in RA; (b) Heat map of DEGs in RA (P was the patient group and C was the control group); (c) Volcano map of DEGs in AMI; (d) Heat map of DEGs in AMI (P was the patient group and C was the control group); (e) Venn diagram of intersection of two samples

图1. (a) RA差异基因火山图；(b) RA差异基因热图(P为疾病组，C为对照组)；(c) AMI差异基因火山图；(d) AMI差异基因热图(P为疾病组，C为对照组)；(e) 两样本交集韦恩图

对AMI数据集DEGs进行GO和KEGG功能富集分析；GO富集分析中与BP相关的主要集中在对脂多糖的反应(response to lipopolysaccharide)、对细菌来源分子的反应(response to molecule of bacterial origin)、细胞因子产生的正向调节(positive regulation of cytokine production)等。与CC相关的主要集中在三级颗粒(tertiary granule)、分泌颗粒膜(secretory granule membrane)、特异性颗粒(specific granule)等。与MF相关的主要集中在模式识别受体活性(pattern recognition receptor activity)、免疫球蛋白结合(immunoglobulin binding)、免疫受体活性(immune receptor activity)等。KEGG主要富集的通路是NF-κB信号通路(NF-kappa

(a) (b) (c) (d)

Figure 2. (a) GO function enrichment analysis of DEGs; (b) KEGG function enrichment analysis of DEGs; (c) GO function enrichment analysis of DEGs; (d) KEGG function enrichment analysis of DEGs

图2. (a) GO富集分析结果；(b) KEGG富集分析结果；(c) GO富集分析结果；(d) KEGG富集分析结果

B signaling pathway)、IL-17信号通路(IL-17 signaling pathway)、TNF信号通路(TNF signaling pathway)等(图2(c)，图2(d))。

3.3. PPI网络模块构建及HubGene鉴定

将交集基因导入STRING数据库，构建综合得分大于0.4的交集基因PPI网络，该网络共由14个节点和6个边组成的PPI网络(图3(a))。筛选出7个基因分别是CCR5、IL1RN、TNFAIP6、CRTAM、RGS1、

(a) (b)

Figure 3. (a) PPI network diagram (The size and color of the circle in the gene represent the importance of the gene); (b) The module with the highest score in Cytoscape’s Degree algorithm

图3. (a) PPI网络(基因中所在圆圈的大小以及颜色深浅代表其基因的重要性)；(b) Cytoscape的Degree算法中得分最高的模块

TREM1、PTX3。通过Cytoscape的Degree算法最终得到CCR5、IL1RN为HubGene (图3(b))。

3.4. HubGene外部数据集验证

筛选出的HubGene在GSE55235、GSE48060数据集验证。以调整后的P < 0.05为筛选条件。结果显示，与对照组相比IL1RN在类风湿关节炎及急性心肌梗死患者中均表达上调(图4(a)， 图4(b))，CCR5在急性心肌梗死患者中表达无显著性差异(图4(c)，图4(d))。

(a) (b) (c) (d)ns为无显著性；^*P < 0.05；^**P < 0.01；^***P < 0.001；^****P < 0.0001。

Figure 4. (a) Expression levels of IL1RN in GSE55235; (b) Expression levels of CCR5 in GSE55235; (c) Expression levels of IL1RN in GSE48060; (d) Expression levels of CCR5 in GSE48060

图4. (a) IL1RN在GSE55235中表达水平；(b) CCR5在GSE55235中表达水平；(c) IL1RN在GSE48060中表达水平；(d) CCR5在GSE48060中表达水平

3.5. 预测信号通路及转录调控因子

以IL1RN表达情况在GSE66360中进行KEGG、GSEA富集分析；KEGG主要富集在破骨细胞转化(Osteoclast differentiation)、NF-κB信号通路(NF-kappa B signaling pathway)、利什曼病(Leishmaniasis)、IL-17信号通路(IL-17 signaling pathway)、吞噬体(Phagosome)等(图5(a))。GSEA主要富集在TNF-α/NF-κB信号通路(TNF-α signaling via NF-κB)、干扰素γ反应(interferon gamma response)、炎症反应(inflammatory response)、补体系统(complement)、移植排斥反应(allograft rejection)等(图5(b))。将HubGene导入TRRUST数据库预测出NFκB1为调控核心基因的转录因子(transcription factor, TF)。NFκB1调控的生物学过程包

(a)(b)

Figure 5. (a) KEGG function enrichment analysis; (b) GSEA function enrichment analysis

图5. (a) KEGG富集分析；(b) GSEA富集分析

括白细胞介素-1介导的信号通路负调控(negative regulation of interleukin-1-mediated signaling pathway)、发热(fever generation)、免疫应答(immune response)、异型细胞间粘附的负调控(negative regulation of heterotypic cell-cell adhesion)、对糖皮质激素的反应(response to glucocorticoid) (表1)。

Table 1. TRRUST predicts transcription factor of IL1RN

表1. TRRUST预测调控核心基因的转录因子

3.6. IL1RN与免疫细胞浸润水平的相关性分析

以IL1RN表达情况高低组使用R语言“GSVA包”进行28种免疫细胞的单样本基因集富集分析(ssGSEA)，结果示IL1RN基因与多种免疫细胞呈正相关(图6)。

ns为无显著性；^*P < 0.05；^**P < 0.01；^***P < 0.001；^****P < 0.0001。

Figure 6. ssGSEA

图6. ssGSEA

3.7. 预测潜在药物

将得到的差异基因输入DGIdb数据库挖掘HubGene潜在相互作用药物，得到与IL1RN相关的药物有甲氨蝶呤(methotrexate)、双醋瑞因(diacerein)、低分子肝素(tinzaparin sodium)、氟哌啶醇(haloperidol decanoate) (表2)。

Table 2. DGIdb predicts potentially interacting drugs

表2. DGIdb数据库预测潜在相互作用药物

4. 讨论

目前对于AMI与RA之间的确切关系仍不明确，随着芯片技术及高通量测序技术的发展，产生出大量基因数据，通过研究基因与相关疾病之间的关系，从而更深层地了解疾病，这已成为当前研究趋势。因此，识别差异表达的HubGene及其分子机制在AMI与RA的发病机制和进展中具有重要作用。

本研究中，通过对AMI和RA数据集中的差异基因分析，GO、KEGG富集分析强调了炎症、免疫途径、趋化因子和细胞因子等在这两种疾病中的重要作用。其中有14个交集基因，分别是IL1RN、CCR5、TNFAIP6、GCSAM、CRTAM、ZNF382、MEST、TREM1、BCL2A1、RGS1、KCNJ15、PTX3、C15orf48、LILRB3；通过对Cytoscape软件构建出交集基因的PPI网络，展现其复杂的互作网络关系，得出IL1RN、CCR5为HubGene。进一步验证IL1RN在验证集中均有差异上调，而CCR5在急性心肌梗死验证集中无差异表达。进一步对IL1RN进行分析发现NF-κB信号通路在两种疾病的发展中起到关键作用。同时在TRRUST数据库预测出转录调控因子NFκB1。ssGSEA分析发现IL1RN与多种免疫细胞高度相关。IL1RN在DGIdb数据库中分析得到甲氨蝶呤、双醋瑞因、低分子肝素、氟哌啶醇四种药物。

IL1RN是编码IL-1Ra的基因 [15] ，IL-1Ra是IL-1家族中的重要组成成分，其家族成员还包括IL-1α、IL-1β [16] 。IL-1具有广泛的生物学效应，如发热、促进免疫应答、参与炎症反应、促进伤口愈合、刺激造血等功能。IL-1ra是特异性IL-1抑制物，其表达数值高低与炎症浸润程度密切相关 [17] 。既往研究发现IL1RN与RA和AMI的发病和病情进展有关，Mukundan Attur等发现携带IL1RN风险单倍型与更严重的RA炎症浸润有关 [18] 。Jinghua Gu等发现IL1RN及其部分相关蛋白在非酒精性脂肪肝诱导RA的过程中起着一定的介导作用 [19] 。Patti G等发现心肌梗死住院患者，IL-1Ra水平与心肌坏死的程度显著相关，并可预测住院事件 [20] 。急性心肌梗死和血流动力学不稳定患者在冠状动脉监护病房期间有心室功能障碍的患者IL-1Ra峰值水平显著高于无严重心室功能障碍的患者 [21] 。

NF-κB信号通路在细胞的炎症反应、免疫应答等过程中起到关键性作用 [22] 。NF-κB的错误调节会引发自身免疫病、慢性炎症以及很多癌症。然而NF-κB信号通路在AMI和RA中的文献报道较少。NFκB1是NF-κB转录因子中的一个成员，虽然NF-κB最常被描述为促进炎症及相关癌症 [23] ，但NFκB1可作为肝细胞癌的肿瘤抑制因子 [24] 。Tyrell Cartwright等发现NFκB1是炎症、衰老和癌症的抑制因子，促进或稳定NFκB1的药物有潜力成为强大的抗炎和抗癌药物 [25] 。然而Sona Margaryan [26] 等发现NFκB1基因甲基化与2型糖尿病(T2DM)患者IL-1Ra水平呈负相关，IL1RN和NFκB1基因启动子的低甲基化可能促进IL-1β/IL-1Ra的产生增加，并调节T2DM的慢性炎症。

总之，本研究基于生物信息学分析了AMI与RA两者之间的相互作用，并对差异基因进行了功能富集分析、PPI网络分析和标志物的筛选。最终发现IL1RN可作为AMI与RA发生发展相互作用的潜在标志物，而IL1RN与NFκB1存在负相关，进而影响NF-κB信号通路。甲氨蝶呤、双醋瑞因、低分子肝素、氟哌啶醇可能作为RA合并AMI的潜在药物。本研究通过基因数据进一步为RA合并AMI的治疗方法开拓了新思路。但仍需要进一步的体内外实验验证与支持。

NOTES

^*第一作者。

^#通讯作者。