Identification of key genes and pathways in atherosclerosis using integrated bioinformatics analysis
BMC Medical Genomics volume 16, Article number: 102 (2023)
Atherosclerosis (AS) is a chronic inflammatory disease that might induce severe cardiovascular events, such as myocardial infarction and cerebral infarction. These risk factors in the pathogenesis of AS remain uncertain and further research is needed. This study aims to explore the potential molecular mechanisms of AS by bioinformatics analyses.
GSE100927 gene expression profiles, including 69 AS samples and 35 healthy controls, were downloaded from Gene Expression Omnibus database and indenfied for key genes and pathways in AS.
A total of 443 differentially expressed genes (DEGs) between control and AS were identified, including 323 down-regulated genes and 120 up-regulated genes. The Gene ontology terms enriched by the up-regulated DEGs were associated with the regulation of leukocyte activation, endocytic vesicle, and cytokine binding, while the down-regulated DEGs were associated with negative regulation of cell growth, extracellular matrix, and G protein-coupled receptor binding. KEGG pathway analysis showed that the up-regulated DEGs were enriched in Osteoclast differentiation and Phagosome, while the down-regulated DEGs were enriched in vascular smooth muscle contraction and cGMP-PKG signaling pathway. Using the modular analysis of Cytoscape, we identified 3 modules mainly involved in Leishmaniasis and Osteoclast differentiation. The GSEA analysis showed the up-regulated gene sets were enriched in the ribosome, ascorbated metabolism, and propanoate metabolism. The LASSO Cox regression analysis showed the top 3 genes were TNF, CX3CR1, and COL1R1. Finally, we found these immune cells were conferred significantly higher infiltrating density in the AS group.
Our data showed the pathway of Osteoclast differentiation and Leishmaniasis was involved in the AS process and we developed a three-gene model base on the prognosis of AS. These findings clarified the gene regulatory network of AS and may provide a novel target for AS therapy.
Atherosclerosis is one of the most common diseases in the cardiovascular system, which threatens human health. The vascular disease causes many vascular complications affecting the cardiovascular system [1, 2]. AS mainly affects the large and medium-sized arteries, and ultimately causes the reduction in arterial inflow or loss of organs, resulting in cells or tissues atrophy/necrosis. As we know, atherosclerotic plaque leads to stenosis and the vulnerability of plaque rupture of a coronary artery, which contributes to the increased risk of cardiovascular events and morbidity and mortality.
The pathological manifestations of AS exhibit the formation of fatty-fibrous plaque, including the adhesion of activated leukocytes, the proliferation of vascular smooth muscle cells (VSMC), and the degradation of the extracellular matrix, angiogenesis, and calcification. Advanced atherosclerotic plaques developed a morphologically fibrous cap composed of VSMCs proliferating and migrating from the vessel wall into the plaque. Our previous studies have reported that pathophysiological factors were involved in the development of AS, including the phenotypic switch of VMSC, the source of oxidative stress, and mitochondrial injury [3, 4].
Atherosclerosis is primarily a chronic inflammatory disease of the arterial wall response to injury, which is mediated by the oxidized low-density lipoprotein (LDL) . In acute artery injury, the endothelial cytokine and the leukocyte adhesion promote the extravasation of monocyte into the sub-endothelial space . Mezentsev reported that desialylation of lipoproteins was an important risk factor in AS and that desialylation of LDL causes the uncontrolled accumulation of lipids by the artery cells . Bezsonov reported that mitochondrial DNA (mtDNA) mutations are associated with AS using leukocytes of atherosclerotic patients . However, the molecular mechanisms involved in the AS process are incompletely understood.
More research will be needed in the public health and medical field and might reduce the incidence and mortality. Given that microarray analysis can provide an enormous amount of data about gene expression, these analysis methods have been used to predict novel prognostic markers and therapeutic targets . In the current study, we downloaded the microarray dataset, GSE100927, from the NCBI-Gene Expression Omnibus database (GEO). The GSE100927 dataset contained 69 AS samples and 35 control samples. Differentially expressed genes (DEGs) were analyzed, and the pathway and Gene ontology enrichment analyses were identified. Furthermore, protein–protein interaction (PPI) networks and modular analysis were constructed to identify hub genes. The module genes were enriched in KEGG pathways. Predictors were identified using the LASSO procedure and inserted into multivariate logistic regression models. We identified the different infiltration of immune cell types between the two groups. These results may provide new insights into the mechanism and pretreatment of AS.
We used the ‘GEOquery’ package of R software (4.0.2 version) to download the GSE100927 microarray data from the GEO database. The GSE100927 data set was based on the GPL17077 platform and contained 104 samples, including 69 atherosclerotic arteries and 35 healthy arteries.
Differentially expressed genes (DEGs) identification
The raw data of GSE100927 was read by the affy package in R (version 3.6.1) and the rate monotonic algorithm was used for the background correction and data normalization. Hierarchical clustering analysis was used to group into 2 similar expression patterns of AS arteries and control. The principal component analysis (PCA) was performed to define probe quality control by the “stats” R package. DEGs were screened by the “limma” package. The cut-off criterion for statistically significant DEGs was determined as |log2FC|> 1 and P value < 0.05. The significant DEGs volcano map was drawn using the “ggplot2” package. The 15 top up-regulated and down-regulated genes were drawn using the “heatmap” package.
Functional enrichment GO and KEGG analysis
Candidate genes enriched functions and pathways were identified by using the online DAVID tool (https://david.ncifcrf.gov/) , which was a free online tool that performed gene annotation and integrated discovery function. Gene ontology analysis (GO) is a common method for the annotation of gene products; it is also able to decide the defining biological characteristics of high-throughput genomic data. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a genetic function database that links functional information [11,12,13]. Both a P < 0.05 and a false discovery rate (FDR) < 0.05 (0.25) were considered to indicate statistical significance.
PPI network and modular analysis
These significant DEGs were conducted in the STRING database (http://string-db.org) to evaluate the relationships of DEGs through protein–protein interaction (PPI) information . The information on protein interaction relationship was imported into the Cytoscape software , which calculated the relationship between DEG encoding proteins in AS. Node degree was analyzed by the Network Analyzer plug-in, which determined the number of interconnections of PPI hub genes. Molecular complex detection (MCODE) was used to select PPI network modules on Cytoscape. The MCODE modules were more than 10 scores which were selected candidates. The nodules of module 1 were conducted for pathway enrichment analyses. Both a P < 0.05 and a false discovery rate (FDR) < 0.05 (0.25) were considered to indicate statistical significance.
Enrichment analysis of gene ontology and KEGG by gene set enrichment analysis (GSEA).
The gene sequences were downloaded from the GSE100927 dataset. GSEA was used to analyze all gene sequences of the samples from two groups , and all genes were sequenced to indicate the trend of gene expression level between the two groups. We conducted GESA analysis and sorted genes according to the algorithm after importing gene annotation files and reference function sets, and then we get a gene sequence table. The GESA software analyzed all gene positions and accumulated them to get enrichment scores.
Screening and verification of diagnostic markers
We used Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression to perform feature selection to screen diagnostic markers for AS arteries . The expression matrix data of GSE100927 were validated as independent matrix data and the diagnostic efficiency was verified as diagnostic markers. The LASSO algorithm was analyzed with the “glmnet” package. The LASSO Cox regression selected the AS features in the model. The LASSO cox regression model was based on the optimal lambda value in the RPI score for each patient. We found the 13 marker genes by using the LASSO cox regression model by the minimum criteria. We combined the genes from LASSO algorithms for further analysis. P < 0.05 was thought to be statistically significant.
The abundance of immune cells in two groups
We analyzed the fractions of immune cells in each sample by CIBERSORT . We performed a Wilcoxon rank-sum test to analyze these immune cells in the two groups. The differential densities of immune cells were shown and the differential colors were represented from low to high infiltrating levels.
Data were expressed as means ± SD. Student’s paired t-test analyses were performed in two groups. The significant differences were considered when P < 0.05. R version 4.0.2 was applied to perform the analysis.
DEGs were separately identified to analyze the microarray data of GSE100927 from each chip. Figure 1A showed the PCA analysis between AS group and the control group. using|logFC|≥ 1 and P < 0.05 as the threshold cutoff points, 443 genes were highlighted. Among these, 323 genes were up-regulated and 120 genes were down-regulated in the AS group compared with the control group. The volcano diagram of 443 DEGs is presented in Fig. 1B. Figure 1C showed a heatmap of the top 15 up-regulated and top 15 down-regulated genes.
GO term enrichment analyses and KEGG pathway analyses
GO and KEGG pathway enrichment analysis of DEGs were performed by using the DAVID online databases. DEGs were classified into biological processes (BP), cellular components (CC), and molecular functions (MF). The top 7 GO enrichment analysis results of DEGs were listed in Fig. 2. The up-regulated genes were mainly enriched in muscle contraction, negative regulation of BMP signaling pathway, and negative regulation of cell growth In BP, while the down-regulated DEGs were significantly enriched in BP such as immune response, inflammatory response, and innate immune response (Fig. 2A). The CC analysis demonstrated that the up-regulated genes were significantly enriched in proteinaceous extracellular matrix, extracellular space, extracellular matrix and region, while the down-regulated genes of DEGs were significantly enriched in MHC class II protein complex, plasma membrane, and endocytic vesicle membrane (Fig. 2B). In addition, MF analysis showed that the up-regulated genes in DEGs were enriched in structural constituent of muscle, fibronectin binding, and heparin binding, while the down-regulated genes were enriched in MHC class II receptor activity, peptide antigen binding, and serine-type endopeptidase activity (Fig. 2C).
Figure 2D showed the significant enrichment pathways of the up-regulated DEGs and down-regulated DEGs analyzed by KEGG analysis. The up-regulated DEGs were significantly enriched in vascular smooth muscle contraction and cGMP-PKG signaling pathway, while the down-regulated DEGs were enriched in Leishmaniasis, Tuberculosis, Phagosome, and Osteoclast differentiation.
Identification of hub genes and modular screening
We employed the STRING online database to construct a protein–protein interaction (PPI) network of DEGs. DEGs were imported into the STRING database and a total of 421 nodes and 3529 edges were obtained with a scoring value > 0.4 (Fig. 3A). The cytoHubba plugin of Cytoscape software was used to identify hub genes in the PPI network. The top 10 hub genes with higher degrees of connectivity were PTPRC, TYROBP, TNF, ITGAM, SPI1, LCP2, CSF1R, TLR2, CCR5, and ITGAX. The higher degree and closeness of the top 10 hub genes were listed in Table 1.
Furthermore, the MCODE plugin of Cytoscape was used to identify gene cluster modules and 6 significant gene cluster modules were obtained. The top 3 key modules from the PPI network were selected: Module 1(significant gene cluster) contained 34 nodes and 361 edges with a score of 21.87. The top 7 genes were ITGAM, TLR1, CCR5, ITGAX, TNF, CD68, and CCR2 (Fig. 3B). Module 2 comprised 41 nodes and 288 edges with a score of 14.40. The top 7 genes were HLA-DPA1, HLA-DQA1, HLA-DBP1, HLA-DQB1, HLA-DRB5, HLA-DQB2, and HLA-DQA2 (Fig. 3C). Module 3 consisted of 40 nodes and 180 edges with a score of 9.23. The top 7 genes were CD53, CYBA, VAMP8, ADAM8, LAIR1, CD36, and STK10 (Fig. 3D).
Identification of KEGG pathway analyses of the top 3 Modules
The top 3 modules with high scores indicated that they may play an important role in the PPI network. KEGG pathway terms of these genes for the top 3 modules were analyzed using DAVID online software. The enrichment analysis results showed that Module 1 enriched 37 KEGG pathways (Additional file 1: Table S1), Module 2 enriched 35 KEGG pathways (Additional file 1: Table S2), and Module 3 enriched 14 KEGG pathways (Additional file 1: Table S3). Module 1 was mainly enriched in cytokine-cytokine receptor interaction, Pertussis, Chemokine signaling pathway, and Tuberculosis (Fig. 4A). Module 2 was mainly enriched in Staphylococcus aureus infection, Antigen processing and presentation, Tuberculosis, and phagosome (Fig. 4B). Module 3 was mainly enriched in phagosome and osteoclast differentiation, respectively (Fig. 4C). Furthermore, we found four common KEGG pathways that were Leishmaniasis, phagosome, Hematopoietic cell lineage, and Osteoclast differentiation in the top 3 modules (Fig. 4D).
Function and enrichment analysis of GSEA.
The DEG expression profile of GSE100927 was uploaded into Gene Set Enrichment Analysis (GSEA) software. GSEA analysis showed that 122/186 gene sets were up-regulated, while 64/186 gene sets were down-regulated between AS group and the control group. Table 2 ranks the most significant enrichments in the up-regulated gene sets (NES > 0.8 and P value < 0.05) and down-regulated gene sets (NES < − 1.3 and P value < 0.05) based on normalized enrichment score (NES). We presented the top 4 most significant plots in the up-regulated (Fig. 5A) and in the down-regulated gene set (Fig. 5B). The up-regulated gene sets were significantly enriched for ribosome, ascorbated and aldarate metabolism, propanoate metabolism, and cardiac muscle contraction, while the down-regulated gene sets were significantly enriched for N-glycanbiosynthesis, rig-I-like receptor signaling pathway, base excision repair, and toll-like receptor signaling pathway, receptively.
Screening and verification of diagnostic markers
To obtain the diagnostic markers genes, Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis was performed based on the gene expression of DEGs. We used the LASSO logistic regression algorithm to screen the diagnostic markers from 443 DEGs. We randomly selected 70% of samples in GSE10927 as the training set construction of the models and the remaining 30% of samples as the test set (internal validation set). A coefficient profile plot was produced against the log (lambda) sequence. Selection of tuning parameter (lambda) in the LASSO model used tenfold cross-validation via minimum criteria. Figure 6A,B show the procedure of screening. We identified the 17 hub genes (Table 3). The 17 hub genes were performed by the multivariate Cox regression analysis to calculate the coefficient of signature, which represented the weight of each gene (Fig. 6C). In order to validate the robustness of the top 3 genes signature, ROC analysis was performed for the validation cohort. The area under the curve (AUC) for TNF, CX3CR1and COL1R1 was 0.96, 0.94, and 0.93, respectively (Fig. 6D). These results indicated that the top 3-genes prognostic model has a high predictive ability of AS.
Differential abundance of infiltrating immune cells in AS
The KEGG pathway enrichment analyses were related to infection, such as Leishmaniasis, and phagosome. To assess the potential associations, we analyzed the different immune cell types between AS group and the control group. The heatmap presented the different immune cell infiltration patterns between the two groups, where the colors displayed the infiltration density (Fig. 7A). In Addition, we utilized the Wilcoxon rank-sum test to compare the difference and discovered these immune cells had a significantly higher infiltrating density in AS groups, including M0 macrophages, T cells gamma delta, memory B cells, as well as activated mast cells (Fig. 7B). On the other hand, we discovered that several immune cells conferred a significantly lower infiltrating density in AS groups, such as monocytes, T cells CD4 memory resting, M2 macrophages, and activated NK cells.
Atherosclerosis (AS), as a metabolic syndrome, posed a serious threat to human health and caused a series of complications affecting the cardiovascular system, including acute myocardial infarction and cerebral infarction . It is well known that AS was initiated by lipid-mediated vascular inflammation, which enhanced continual monocyte recruitment and foam-cell formation in the arterial wall. Roy reported that AS is characterized by imbalanced lipid metabolism . The GSE100927 data set successfully extracted 443 DEGs between AS and the control group. Functional annotations of DEGs are mainly assessed through the analysis of GO and KEGG pathway enrichment.
BP analyses demonstrated that the up-regulated DEGs were mainly enriched in muscle contraction and negative regulation of the BMP signaling pathway, while the down-regulated DEGs were enriched in immune response and inflammatory response . These findings are in line with the hypothesis, which stated that AS development is a complicated inflammatory process with immune reactions and is associated with decreased vascular smooth muscle contraction and proliferation [5, 6]. In addition, the enriched KEGG pathways of up-regulated DEGs found osteoclast differentiation and phagosome. Previous studies have shown that lipid phagocytosis is a crucial process involved in atherogenesis [19, 20] and the inhibition of osteoclastic differentiation impaired calcium resorption in atherosclerotic plaques . Our result was consistent with these findings that the up-regulated genes of DEGs were enriched in vascular smooth muscle contraction, while the down-regulated genes were enriched in Leishmaniasis, Phagosome, and Osteoclast differentiation.
We identified the top 10 hub genes based on degree value by using Cytoscape software, including PTPRC, TYROBP, TNF, ITGAM, SPI1, LCP2, CSF1R, CCR5, and ITGAX. Nie reported that PTPRC was considered to be an immune marker gene in atherosclerosis development . Liu reported that the upregulated TYROBP was identified as a common hub gene and as an immune-responsive hub gene in the advanced atherosclerotic plaques . Lu found that TYROBP expression was significantly higher in the low-grade glioma tissues compared to the normal tissues . TNF-αis a pro-atherogenic stimulus on vascular endothelial cells. Yuan supported the evidence of the causal associations of increased TNF levels with a higher risk of cardiovascular diseases, including AS . The TNF-α expression was higher in obese patients with cardiovascular disease (CAD) . Wei reported that Csf1r inhibited macrophage proliferation and impaired the progression of AS . Lin reported that CCL5/CCR5 axis activation promoted the proliferation and the phenotypic switching of vascular smooth muscle cells . A study found CCR5 is related to better survival of triple-negative breast cancer through improvement tumor microenvironment . Zhao reported that ITGAM is a potential diagnostic and prognostic biomarker of unstable atherosclerotic plaque-related stroke, which will be a novel therapeutic target and the development of novel management strategies . Therefore, these hub genes are closely related to the pathological process of AS and may play a key role in the early diagnosis of AS patients. The PPI network construction highlighted these important genes, which revealed the characteristics of atherosclerotic lesions.
The functions of the genes in three important modules from PPI were selected to further perform the KEGG pathway enrichment analysis. The common three KEGG enrichment pathways were osteoclast differentiation, phagosome, and Leishmaniasis. These findings are in line with other research which reported that AS development is a result of imbalanced lipid metabolism and the AS progression was mediated by the phagosome of macrophages and smooth muscle cells . The findings were consistent with our previous research which reported that the transdifferentiation of vascular smooth muscle cells mediated the AS process . The GSEA results showed that these genes were enriched in ribosome, N-glycanbiosynthesis, and rig-I-like receptor signaling pathway, which were in line with the KEGG pathway analysis identified phagosome, Leishmaniasis, and osteoclast differentiation. Some hub genes, including TYROBP, TNF, SP11, LCP2 and CSF1R, played important roles in phagosome [31, 32]. TNF and SP11 are related to leishmaniasis infection [33, 34]. Moreover, CCR5 is strongly associated with Leishmaniasis in HIV-1-infected patients  and the CCR5 reduction accelerated osteoclastogenesis . Some studies found that these hub genes, such as PTPRC, TNF, ITGAM, SP11, and ITGAX, are significantly associated with osteoclast differentiation [34, 37, 38]. CSF1R expression elucidated the impact of specific miRNA-mRNA interactions in osteoclastogenesis , All in all, these hub genes are highly connected to three enriched pathways.
To evaluate the prognostic genes based on DEGs of GSE100927, we applied the LASSO Cox regression model to construct a prognostic gene signature. We identified a gene signature (TNF, CX3CR1, and COL1R1) involved in AS. The AUC value of the three-gene signature suggested the prediction ability and the prognostic values of these gene signatures. TNF-α could induce the up-expression of CX3CR1 in endothelial cells and BV-2 microglial cells [40, 41]. Morimura reported that the peritoneal macrophages of the CX3CR1 (−/−) mice expressed lower expression levels of TNF-α, IL-1β, and IL-6 . Garre suggested the highlight CX3CR1high monocytes and TNF-α as potential therapeutic targets for preventing infection-induced cognitive dysfunction . The novel insight provides new noninvasive methods and may inform early prediction and early intervention strategies for AS patients.
AS is a chronic inflammatory disease involving multiple types of immune cells. Chistiakov reported that Tregs mediated the immune response and the secretion of anti-inflammatory cytokines IL-10 and TNF-beta . CD8 + T cells enhanced the secretion of IFN-γ, IL-2, and TNF-α . DC promoted the activating of the TNF-α/NF-κB/CXCR-4 pathway . Ni reported that CX3CR1 deficiency caused a significant increase in inflammatory monocyte/macrophage infiltration . Collagen type I(COL1) is remodeled in atherosclerotic plaque and the COL1 degradation is associated with mortality in patients with AS . These results were consistent with the immune cell infiltration in AS groups.
The main pathological process of AS involved vascular endothelial cell damage, lipid deposition, and plaque formation . Our results showed that the down-regulated genes enriched in Leishmaniasis, Phagosome, and Osteoclast differentiation. We identified the top 10 hub genes. Some hub genes(TYROBP, TNF, SP11, LCP2, and CSF1R) played a key role in phagosome. TNF, SP11, and CCR5 are related to leishmaniasis infection. Other hub genes(PTPRC, TNF, ITGAM, SP11, and ITGAX) are significantly associated with osteoclast differentiation. We constructed a gene signature (TNF, CX3CR1, and COL1R1). The AUC value suggested the prediction ability and the prognostic value. Three genes were associated with immune cell infiltration in AS patients. Our findings explored the regulatory network of AS and represented a novel diagnostic and therapeutic target for AS.
There are also some limitations in our study. Firstly, the hub genes and the gene signature were analyzed using gene expression profiles of human samples from online public databases; further studies using more human samples will be required to validate these marker genes by WB experiments. Secondly, more experiments need to examine these potential mechanisms related to the phagosome, Leishmaniasis, and osteoclast differentiation in tissue samples from patients with AS.
In summary, our study displays a new comprehensive bioinformatics analysis of DEGs that might be involved in the AS progress. We identified 10 hub genes and three modules that were highly associated with AS. The KEGG pathways were mainly enriched in phagosome, Leishmaniasis, and Osteoclast differentiation. We developed and validated a three-gene signature to predict the prognosis of AS. These novel molecular targets are principal candidate genes for further investigations into biomarkers and molecular mechanisms.
Availability of data and materials
The datasets analyzed during the current study are available in the GEO repository, (https://www.ncbi.nlm.nih.gov/geo, accession number GSE100927)..
Gene Expression Omnibus
Kyoto Encyclopedia of Genes and Genomes
Differentially expressed genes
Least absolute shrinkage and selection operator
Vascular smooth muscle cell
Principal component analysis
False discovery rate
Molecular complex detection
Search Tool for the Retrieval of Interacting Genes/Proteins
Gene set enrichment analysis
Major histocompatibility complex
Tumor necrosis factor-α
CX3C receptor 1
Lanzer P, Hannan F, Lanzer J, Janzen J, Raggi P, Furniss D, et al. Medial arterial calcification: JACC state-of-the-art review. J Am Coll Cardiol. 2021;78(11):1145–65.
Xu S, Ilyas I, Little P, Li H, Kamato D, Zheng X, et al. Endothelial dysfunction in atherosclerotic cardiovascular diseases and beyond: from mechanism to pharmacotherapies. Pharmacol Rev. 2021;73(3):924–67.
Durham A, Speer M, Scatena M, Giachelli C, Shanahan C. Role of smooth muscle cells in vascular calcification: implications in atherosclerosis and arterial stiffness. Cardiovasc Res. 2018;114(4):590–600.
Wesseling M, Sakkers T, de Jager S, Pasterkamp G, Goumans M. The morphological and molecular mechanisms of epithelial/endothelial-to-mesenchymal transition and its involvement in atherosclerosis. Vascul Pharmacol. 2018;106:1–8.
Roy P, Orecchioni M, Ley K. How the immune system shapes atherosclerosis: roles of innate and adaptive immunity. Nat Rev Immunol. 2021. https://doi.org/10.1038/s41577-021-00584-1.
Libby P. The changing landscape of atherosclerosis. Nature. 2021;592(7855):524–33.
Mezentsev A, Bezsonov E, Kashirskikh D, Baig MS, Eid AH, Orekhov A. Proatherogenic sialidases and desialylated lipoproteins: 35 years of research and current state from bench to bedside. Biomedicines. 2021. https://doi.org/10.3390/biomedicines9060600.
Bezsonov EE, Sobenin IA, Orekhov AN. Immunopathology of atherosclerosis and related diseases: focus on molecular biology. Int J Mol Sci. 2021. https://doi.org/10.3390/ijms22084080.
Chen J, Zhang X, Millican R, Sherwood J, Martin S, Jo H, et al. Recent advances in nanomaterials for therapy and diagnosis for atherosclerosis. Adv Drug Deliv Rev. 2021;170:142–99.
da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Kanehisa M, Sato Y, Kawashima M. KEGG mapping tools for uncovering hidden features in biological data. Protein Sci. 2022;31(1):47–53.
Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51.
Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021;49(D1):D605-d612.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Subramanian A, Tamayo P, Mootha V, Mukherjee S, Ebert B, Gillette M, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102(43):15545–50.
Efendi A, Ramadhan HW, editors. Parameter estimation of multinomial logistic regression model using least absolute shrinkage and selection operator (LASSO). In: The 8th annual basic science international conference: coverage of basic sciences toward the World’s Sustainability Challanges (2018).
Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37(7):773–82.
Zhao TX, Mallat Z. Targeting the immune system in atherosclerosis: JACC state-of-the-art review. J Am Coll Cardiol. 2019;73(13):1691–706.
Moore KJ, Koplev S, Fisher EA, Tabas I, Bjorkegren JLM, Doran AC, et al. Macrophage trafficking, inflammatory resolution, and genomics in atherosclerosis: JACC macrophage in CVD series (Part 2). J Am Coll Cardiol. 2018;72(18):2181–97.
Zavaczki E, Gall T, Zarjou A, Hendrik Z, Potor L, Toth CZ, et al. Ferryl hemoglobin inhibits osteoclastic differentiation of macrophages in hemorrhaged atherosclerotic plaques. Oxid Med Cell Longev. 2020;2020:3721383.
Nie H, Yan C, Zhou W, Li T. Analysis of immune and inflammation characteristics of atherosclerosis from different sample sources. Oxid Med Cell Longev. 2022;2022:5491038.
Liu C, Zhang H, Chen Y, Wang S, Chen Z, Liu Z, et al. Identifying RBM47, HCK, CD53, TYROBP, and HAVCR2 as hub genes in advanced atherosclerotic plaques by network-based analysis and validation. Front Genet. 2020;11: 602908.
Lu J, Peng Y, Huang R, Feng Z, Fan Y, Wang H, et al. Elevated TYROBP expression predicts poor prognosis and high tumor immune infiltration in patients with low-grade glioma. BMC Cancer. 2021;21(1):723.
Yuan S, Carter P, Bruzelius M, Vithayathil M, Kar S, Mason AM, et al. Effects of tumour necrosis factor on cardiovascular disease and cancer: a two-sample Mendelian randomization study. EBioMedicine. 2020;59: 102956.
Bilgic Gazioglu S, Akan G, Atalar F, Erten G. PAI-1 and TNF-alpha profiles of adipose tissue in obese cardiovascular disease patients. Int J Clin Exp Pathol. 2015;8(12):15919–25.
Wei Y, Zhu M, Corbalan-Campos J, Heyll K, Weber C, Schober A. Regulation of Csf1r and Bcl6 in macrophages mediates the stage-specific effects of microRNA-155 on atherosclerosis. Arterioscler Thromb Vasc Biol. 2015;35(4):796–803.
Lin CS, Hsieh PS, Hwang LL, Lee YH, Tsai SH, Tu YC, et al. The CCL5/CCR5 axis promotes vascular smooth muscle cell proliferation and atherogenic phenotype switching. Cell Physiol Biochem. 2018;47(2):707–20.
Wang X, Han Y, Peng J, He J. CCR5 is a prognostic biomarker and an immune regulator for triple negative breast cancer. Aging (Albany NY). 2021;13(20):23810–30.
Zhao B, Wang D, Liu Y, Zhang X, Wan Z, Wang J, et al. Six-gene signature associated with immune cells in the progression of atherosclerosis discovered by comprehensive bioinformatics analyses. Cardiovasc Ther. 2020;2020:1230513.
Xu X, Hao Y, Wu J, Zhao J, Xiong S. Assessment of weighted gene co-expression network analysis to explore key pathways and novel biomarkers in muscular dystrophy. Pharmgenom Pers Med. 2021;14:431–44.
Brawn LC, Hayward RD, Koronakis V. Salmonella SPI1 effector SipA persists after entry and cooperates with a SPI2 effector to regulate phagosome maturation and intracellular replication. Cell Host Microbe. 2007;1(1):63–75.
Salih MAM, Fakiola M, Lyons PA, Younis BM, Musa AM, Elhassan AM, et al. Expression profiling of Sudanese visceral leishmaniasis patients pre-and post-treatment with sodium stibogluconate. Parasite Immunol. 2017. https://doi.org/10.1111/pim.12431.
Yao Z, Getting SJ, Locke IC. Regulation of TNF-induced osteoclast differentiation. Cells. 2021. https://doi.org/10.3390/cells11010132.
Vallejo A, Abad-Fernandez M, Moreno S, Moreno A, Perez-Elias MJ, Dronda F, et al. High levels of CD4(+) CTLA-4(+) Treg cells and CCR5 density in HIV-1-infected patients with visceral leishmaniasis. Eur J Clin Microbiol Infect Dis. 2015;34(2):267–75.
Lee D, Shin KJ, Kim DW, Yoon KA, Choi YJ, Lee BNR, et al. CCL4 enhances preosteoclast migration and its receptor CCR5 downregulation by RANKL promotes osteoclastogenesis. Cell Death Dis. 2018;9(5):495.
Li Z, Sun Y, He M, Liu J. Differentially-expressed mRNAs, microRNAs and long noncoding RNAs in intervertebral disc degeneration identified by RNA-sequencing. Bioengineered. 2021;12(1):1026–39.
Skinkyte-Juskiene R, Kogelman LJA, Kadarmideen HN. Transcription factor co-expression networks of adipose RNA-Seq data reveal regulatory mechanisms of obesity. Curr Genom. 2018;19(4):289–99.
Hrdlicka HC, Lee SK, Delany AM. MicroRNAs are critical regulators of osteoclast differentiation. Curr Mol Biol Rep. 2019;5(1):65–74.
Szukiewicz D, Wojciechowska M, Bilska A, Stangret A, Szewczyk G, Mittal TK, et al. Aspirin action in endothelial cells: different patterns of response between chemokine CX3CL1/CX3CR1 and TNF-alpha/TNFR1 signaling pathways. Cardiovasc Drugs Ther. 2015;29(3):219–29.
Li X, Ye Z, Guo Q, Wang E, Pan Y. PDTC ameliorates neuropathic pain by inhibiting microglial activation via blockage of the TNFalpha-CX3CR1 pathway. Eur J Histochem. 2021. https://doi.org/10.4081/ejh.2021.3184.
Morimura S, Oka T, Sugaya M, Sato S. CX3CR1 deficiency attenuates imiquimod-induced psoriasis-like skin inflammation with decreased M1 macrophages. J Dermatol Sci. 2016;82(3):175–88.
Garre JM, Silva HM, Lafaille JJ, Yang G. CX3CR1(+) monocytes modulate learning and learning-dependent dendritic spine remodeling via TNF-alpha. Nat Med. 2017;23(6):714–22.
Chistiakov DA, Sobenin IA, Orekhov AN. Regulatory T cells in atherosclerosis and strategies to induce the endogenous atheroprotective immune response. Immunol Lett. 2013;151(1–2):10–22.
Yang J, Liu R, Deng Y, Qian J, Lu Z, Wang Y, et al. MiR-15a/16 deficiency enhances anti-tumor immunity of glioma-infiltrating CD8+ T cells through targeting mTOR. Int J Cancer. 2017;141(10):2082–92.
Han N, Li X, Wang Y, Wang L, Zhang C, Zhang Z, et al. Increased tumor-infiltrating plasmacytoid dendritic cells promote cancer cell proliferation and invasion via TNF-alpha/NF-kappaB/CXCR-4 pathway in oral squamous cell carcinoma. J Cancer. 2021;12(10):3045–56.
Ni Y, Zhuge F, Ni L, Nagata N, Yamashita T, Mukaida N, et al. CX3CL1/CX3CR1 interaction protects against lipotoxicity-induced nonalcoholic steatohepatitis by regulating macrophage migration and M1/M2 status. Metabolism. 2022;136: 155272.
Li B, Song X, Guo W, Hou Y, Hu H, Ge W, et al. Single-cell transcriptome profiles reveal fibrocytes as potential targets of cell therapies for abdominal aortic aneurysm. Front Cardiovasc Med. 2021;8: 753711.
This work was supported by grants from the Fund of the National Natural Science Foundation of China (No. 82073852) and the Science and Technology Department of Hubei Province (ZX2021F007) for the design of the study and the Found of Hubei University of Science and Technology (2018xzy02, 2019xz01, 2019-21GP12) for the analysis and in writing the manuscript. No specific funding was received for this study.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The KEGG pathway of Module 1, Related to Figure 4A. Table S2. The KEGG pathway of Module 2, Related to Figure 4B. Table S3. The KEGG pathway of Module 3, Related to Figure 4C.
About this article
Cite this article
Li, S., Li, S., Li, Q. et al. Identification of key genes and pathways in atherosclerosis using integrated bioinformatics analysis. BMC Med Genomics 16, 102 (2023). https://doi.org/10.1186/s12920-023-01533-8
- Differentially expressed genes (DEGs)
- Network modules
- Bioinformatics analysis