Identification of key genes and biological pathways associated with vascular aging in diabetes based on bioinformatics and machine learning

Data download and preprocessing

Four gene expression datasets (GSE66280 and GSE171663, GSE121487 and GSE57329) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geoprofiles/). GSE66280 and GSE171663 are gene expression datasets related to vascular smooth muscle cell (VSMC) senescence and were used to identify genes associated with VSMC senescence. GSE66280 is based on the GPL16570 platform (Affymetrix Mouse Gene 2.0 ST Array) and includes six normal VSMC samples and six VSMC senescence samples induced by high glucose. GSE171663 is based on the GPL18573 platform (Illumina NextSeq 500) and includes four normal VSMC samples and four stress-induced VSMC senescence samples. GSE121487 and GSE57329 targeted diabetic vasculopathy, both incorporating GPL7202 platform comprising gene expression from three normal and three diabetic mouse aorta samples each.

The limma R package was used to perform background correction, median normalization, and gene symbol conversion on the four datasets. The GSE121487 and GSE57329 datasets were merged and batch effects were removed using the “sva” package. Gene expression box plots and Uniform Manifold Approximation and Projection (UMAP) analysis were performed to verify the effectiveness of batch correction.

VSMC senescence-related genes screening

The R package “limma” was used to perform differential expression analysis between normal VSMC samples (Control group) and senescent VSMC samples (Senescence group) in the GSE66280 and GSE171663 datasets, respectively. Genes with a threshold of “|fold change| >1.2 and p < 0.05” were selected as differentially expressed genes (DEGs). The common DEGs in the two datasets were defined as VSMC senescence related genes.

Diabetic vascular aging-related genes screening

Differential analysis of gene expression was performed between normal mouse aorta samples (Control group) and diabetic mouse aorta samples (Diabetes group) in the merged GSE121487 and GSE57329 datasets related to diabetic vascular aging using the R package “limma”. DEGs were selected based on the threshold of “|fold change| >1.2 and p < 0.05”. By comparing the genes related to VSMC senescence, the common DEGs were identified as genes that affect diabetic vascular aging.

Functional enrichment analysis and protein-protein interaction (PPI) analysis

Functional enrichment analysis, including KEGG (Kyoto Encyclopedia of Genes and Genomes) and Gene Ontology (GO) annotation, is a bioinformatics approach that associates specific sets of genes with known signaling pathways and biological functions to identify the signaling pathways and biological functions involved in the target genes. The GO annotation includes three categories: biological processes (BPs), cellular components (CCs), and molecular functions (MFs). In this study, we used the DAVID database (https://david.ncifcrf.gov) and Metascape dataset (http://metascape.org/gp/index.html#/main/step1) to perform KEGG and GO functional enrichment analysis for the VSMC senescence-related genes and diabetic vascular aging-related genes, respectively. For DAVID analysis, the functional annotation tool was used with the following settings: Species: Mus musculus, Annotation Categories: KEGG_PATHWAY, GO_TERM_BP_DIRECT, GO_TERM_CC_DIRECT, GO_TERM_MF_DIRECT, and the default settings for remaining parameters. For Metascape analysis, the “Express Analysis” tool was used with the following parameters: Species: Mus musculus, Analysis Mode: Custom Analysis, and the default settings for remaining parameters. A p-value < 0.05 was considered statistically significant in these functional enrichment analyses. In addition, to identify significantly different signaling pathways between two different biological conditions in the dataset, we chose KEGG gene sets genomes as the reference genome and performed Gene Set Enrichment Analysis (GSEA) using GSEA software version 4.0.1 and gene set variation analysis (GSVA) using the “GSVA” R package. A p-value < 0.05 was considered significant in these functional enrichment analysis results.

To explore the interactions among the diabetic vascular aging-related genes, we uploaded these genes to the STRING database (https://string-db.org/). A PPI network model was constructed using Cytoscape 3.7.1 software, and the topological parameters of each target were calculated.

Selection of key genes based on machine learning

Two machine learning algorithms, two robust machine learning algorithms, Least Absolute Shrinkage and Selection Operator (LASSO) regression and Support Vector Machine-Recursive Feature Elimination (SVM-RFE), were used to select key genes related to diabetic vascular aging from the diabetic vascular aging dataset. LASSO regression, implemented using the “glmnet” R package, was first applied to the set of genes related to diabetic vascular aging. The optimal regularization strength (λ) was determined through 10-fold cross-validation, and the LASSO model was trained on the entire dataset using the selected λ value. Genes with non-zero coefficients above a threshold of |coefficient| >0.001 were considered important features. In parallel, SVM-RFE analysis was conducted using the “e1071” R package, which iteratively eliminates the least important features from the SVM model based on the absolute values of the SVM weight vector. A grid search combined with 5-fold cross-validation was performed to determine the optimal number of features to retain. The intersection genes obtained from the two analyses were considered as key genes related to vascular aging in diabetes.

Immune infiltration analysis

CIBERSORTx (https://cibersortx.stanford.edu/) is a computational tool based on single-cell RNA sequencing (scRNA-seq) data that can be used to infer the types and relative proportions of immune cells in tissue samples. In this study, the gene expression data of each sample in the diabetic vascular aging dataset was uploaded to the CIBERSORTx database to calculate the infiltration levels of 22 immune cell types. Wilcoxon rank-sum test was used to compare the differences in immune cell infiltration between normal vascular tissue and diabetes vascular tissue, and a p-value < 0.05 was considered significant. Correlations were then examined between immune cells, as well as between key genes and immune cells.

Expression of key genes and receiver operating characteristic (ROC) analysis

The Wilcoxon rank-sum test was used to analyze the expression levels of key genes in GSE66280, GSE171663, and the merged diabetic vasculopathy dataset (merged GSE121487 and GSE57329 datasets). To evaluate the specificity and sensitivity of these key genes in the diagnosis of diabetic vascular aging, we conducted receiver operating characteristic (ROC) analysis on GSE66280, GSE171663, and the merged diabetic vasculopathy dataset, and determined the area under the curve (AUC) values using the “pROC” package.

Availability of data and materials

The entire gene expression profile data in this study were downloaded from the Gene Expression Omnibus (GEO) database (Website: https://www.ncbi.nlm.nih.gov/geo/).

Screening and functional enrichment analysis of VSMC senescence-related genes

A total of 4085 DEGs were identified between senescent VSMC samples (Senescence group) and normal VSMC samples (Control group) in the GSE66280 dataset, including 1882 upregulated genes and 2203 downregulated genes in the Senescence group (Figure 1A). In the GSE171663 dataset, a total of 5110 DEGs were identified, including 1981 upregulated genes and 3129 downregulated genes (Figure 1B). Through a comparison of the two datasets, 428 common DEGs were identified, including 126 upregulated genes and 302 downregulated genes in the Senescence group (Figure 1C). These common DEGs were considered to be related to VSMC senescence. The expression of these VSMC senescence related genes in the GSE66280 and GSE171663 datasets was shown in the heat map (Figure 1D, 1E).

Functional enrichment analysis was performed on these common DEGs. GO annotation results showed that these genes exerted multiple functions, such as protein binding, RNA binding, ATP binding, rRNA binding, and pseudouridine synthase activity. These genes were mainly enriched in many cellular components (CCs), such as cytosol, nucleus, nucleoplasm, cytoplasm, membrane. And they were involved in several BPs such as rRNA processing, ribosomal large subunit biogenesis, mitochondrial translation, ribosome biogenesis, and regulation of cell cycle. The top 10 most significantly enriched GO terms were displayed in Figure 2A. KEGG pathway enrichment analysis showed that vascular smooth muscle cell senescence-related genes were enriched in 17 signaling pathways, including ribosome biogenesis in eukaryotes, ribosome, cellular senescence, biosynthesis of cofactors, PI3K-Akt signaling pathway (Figure 2B, 2C). We further analyzed the expression differences of these KEGG pathways between different groups in the GSE66280 and GSE171663 datasets using gene set enrichment analysis (GSEA). The results showed that in the GSE66280 dataset, 9 KEGG pathways were significantly different between normal and senescent VSMCs. Six KEGG pathways, such as cellular senescence, ECM-receptor interaction, PI3K-Akt signaling pathway, and the AGE-RAGE signaling pathway in diabetic complications, were activated in senescent VSMCs, while three KEGG pathways, including ribosome biogenesis in eukaryotes, biosynthesis of cofactors, and ribosome, were inhibited in senescent VSMCs (Figure 2D). In the GSE171663 dataset, 4 KEGG pathways, including cellular senescence, ECM-receptor interaction, focal adhesion, and apoptosis, were activated in senescent VSMCs, while the ribosome pathway was inhibited (Figure 2E).

Screening and functional enrichment analysis of diabetic vascular aging-related genes

The GSE121487 and GSE57329 expression matrices were combined and batch effects were removed. Gene expression boxplot and principal component analysis were performed on the corrected data. As shown in the gene expression boxplot, the data distribution between the different datasets became consistent after batch effect removal, with the medians aligned on a single line (Figure 3A). UMAP analysis showed that samples from GSE121487 and GSE57329 were randomly distributed, indicating that the influence of batch effect had been eliminated (Figure 3B). After batch effect removal, a total of 1604 differentially expressed genes were identified between diabetes vascular samples (Diabetes group) and normal vascular samples (Control group) (Figure 3C). By comparing with the screened VSMC senescence-related genes, 52 differentially expressed genes were obtained, including 16 significantly upregulated genes in the Diabetes group and 36 significantly downregulated genes (Figure 3D, Supplementary Table 1). These genes were considered as diabetic vascular aging-related genes. The expression of these diabetic vascular aging-related genes in the diabetes group and control group was shown in the heat map (Figure 3E).

The interactions among these diabetic vascular aging-related genes were visualized in a PPI network. Among them, CCT2, RPS11, RPS3, NPM1, EIF3G, TFB1M and DALRD3 were identified as the genes with the highest “degree” values in the PPI network, suggesting their pivotal roles in the regulation of diabetic vascular aging-related genes (Figure 3F). Functional enrichment analysis of these genes was performed using the Metascape database. The results revealed that these genes were mainly enriched in 41 GO terms and 8 KEGG signaling pathways, such as ribosomal subunit, AMPK signaling pathway, unfolded protein binding, ribonucleoprotein complex biogenesis, oxidoreductase activity, and RNA modification. To visualize the interactions, the GO terms and KEGG pathways with the strongest interconnections were selected for network analysis (Figure 3G). Furthermore, the GSVA results showed that 30 KEGG signaling pathways, such as the beta-alanine metabolism, ribosome, mammalian circadian rhythm, endocytosis and cell adhesion molecules (CAMs), were significantly upregulated in diabetic vascular aging. In contrast, 10 KEGG signaling pathways, such as fatty acid metabolism, VEGF signaling, PPAR signaling and peroxisome, were significantly downregulated in diabetes-associated vascular aging (Figure 3H).

Selection of key genes based on machine learning

Two machine learning algorithms were used in this study to further screen for key genes associated with vascular aging in diabetes. These key genes could serve as potential biomarkers or therapeutic targets for diabetes associated vascular aging. First, a LASSO regression model was designed based on diabetes vascular samples and normal vascular samples. When λ was set to 0.04, the LASSO regression model accurately distinguished between diabetes vascular samples and normal vascular samples, yielding 12 candidate genes (Figure 4A). Moreover, SVM–REF analysis was performed and identified 15 feature genes (Figure 4B). Finally, the results from the two algorithms were combined, yielding 6 key genes associated with diabetic vascular aging, including TFB1M, FOXRED2, LY75, DALRD3, PI4K2B, and NDOR1 (Figure 4C).

Immune infiltration analysis

We analyzed the immune cell infiltration in diabetic vascular tissues and normal vascular tissues. A total of 19 infiltrating immune cell types were identified in the vascular tissues. The bar graph displayed the type and number of immune cell infiltrates in each sample (Figure 5A). Subsequently, we assessed the correlation between these immune cell populations. A strong positive correlation between Mast cells activated and B cells memory (r = 0.9) and a strong negative correlation between T cells CD4 naive and Eosinophils (r = −0.78) was found (Figure 5B).

In addition, we compared the differences in infiltration levels of the 19 immune cell types between diabetic vascular tissues and normal vascular tissues. Compared to normal vascular tissues, T cells CD4 naive and Macrophages M1 showed significantly higher expression in diabetic vascular tissues, whereas T cells CD8, Eosinophils, and T cells regulatory (Tregs) exhibited significantly lower expression in diabetic vascular tissues (Figure 5C).

Furthermore, we investigated the correlation between the expression of six key genes (TFB1M, FOXRED2, DALRD3, PI4K2B, and NDOR1) and the infiltration of the 19 immune cell types (Figure 6A–6F). The differentially expressed immune cells in diabetic vascular tissues and normal vascular tissues showed a certain correlation with the expression of these key genes. Specifically, TFB1M expression displayed a significant positive association with the infiltration status of T cells CD4 naive and negatively correlated with the infiltration status of Eosinophils and T cells regulatory (Tregs). PI4K2B expression showed a significant positive link to the infiltration status of T cells CD4 naive and Macrophages M1 and negatively correlated with T cells CD8 and Eosinophils. FOXRED2 expression was significantly positively linked to the infiltration status of T cells CD4 naive and Macrophages M1 and negative correlation with the infiltration status of T cells CD8. NDOR1 expression exhibited a significantly positive correlation with naive CD4 T cell infiltration levels and a significantly negative correlation with Eosinophil infiltration levels. DALRD3 expression showed a significant positive correlation with T cells CD4 naive infiltration levels. No significant correlation was found between LY75 and immune cell infiltration.

Expression of key genes and ROC analysis

In datasets GSE66280 (Figure 7A) and GSE171663 (Figure 7B), the expression of key genes TFB1M, FOXRED2, DALRD3, PI4K2B, and NDOR1 was significantly downregulated in senescent VSMCs compared to normal VSMCs, while the expression of LY75 gene was significantly upregulated. Similarly, in the merged diabetic vasculopathy dataset (GSE121487 and GSE57329), TFB1M, FOXRED2, DALRD3, PI4K2B and NDOR1 were markedly downregulated, whereas LY75 was upregulated in diabetic blood vessels (Figure 7C).

The ROC analysis revealed that these key genes can accurately distinguish between senescent VSMCs and normal VSMCs, as well as between normal vascular samples and diabetic vascular aging samples. The AUC values of TFB1M, FOXRED2, LY75, DALRD3, PI4K2B, and NDOR1 were all greater than 0.8 across all analyzed datasets, indicating their strong potential as biomarkers for vascular aging (Figure 8A–8C).