Weighted gene co-expression network analysis (WGCNA) and module screening
GLN-associated gene sets were investigated using WGCNA. The results demonstrated that when the weighted value was 24 (Fig. 1A), scale independence was greater than 0.85 and mean connectivity was approximately 0. Three co-expressed modules were screened and the unrelated genes were distributed to a gray module, which was excluded from subsequent analyses (Fig. 1B). The module eigengenes (MEs) were correlated, thus allowing studying the associations among modules. A dendrogram and heatmap were adopted to depict the eigengene network (Fig. 1C). To understand the physiological significance of these modules, three MEs were associated with GLN, and the most significant associations were searched based on the heatmap of module-trait correlations (Fig. 1D).
A total of 1039 statistically significant DEGs were identified between healthy controls and atherosclerosis samples (p-adjusted < 0.05, |log2 fold change [FC]|> 0.5). Atherosclerosis samples included 619 genes with upregulated and 420 with downregulated expression. A volcano plot was generated for DEG visualization (Fig. 2A). The top 5 upregulated (WDFY4, TPP1, DAB2, CYTH4, and ADAP2) and downregulated DEGs (ATP1A2, CAB39L, BAG2, C14orf132, and PXYLP1) are shown using a heatmap (Fig. 2B). Based on the results of the Wilcoxon tests, these 10 genes showed marked differences in expression levels between atherosclerosis and control samples (p < 0.05, Fig. 2C).
We obtained 308 GLN-associated DEGs from the intersection between GLN-associated module genes and DEGs.
To explore potential mechanisms underlying DEG function, GSEA was performed. A Molecular Signatures Database (MSigDB) collection was used to select the signaling pathways with the most significant enrichment according to the normalized enrichment score (NES). GSEA revealed that LYSOSOME (NES = 2.45, P-adjusted = 0.012, FDR = 0.008), HEMATOPOIETIC CELL LINEAGE (NES = 2.41, P-adjusted = 0.012, FDR = 0.008), CYTOKINE–CYTOKINE RECEPTOR INTERACTION (NES = 2.351, P-adjusted = 0.012, FDR = 0.008), PROPANOATE METABOLISM (NES = − 1.795, P-adjusted = 0.023, FDR = 0.014), DILATED CARDIOMYOPATHY (NES = − 1.796, P-adjusted = 0.015, FDR = 0.009), and TYROSINE METABOLISM (NES = − 1.798, P-adjusted = 0.014, FDR = 0.009) (Fig. 3A–F) were significantly enriched in atherosclerosis.
Gene set variation analysis (GSVA)
To further explore functional annotation, differences in the relative expression of pathways were evaluated between disease and control groups using GSVA. Numerous differentially expressed pathways were enriched and a heatmap was used for visualization. The disease group exhibited markedly lower expression of KEGG_LYSOSOME- and KEGG_HEMATOPOIETIC_CELL_LINEAGE-associated pathways, and markedly higher expression of KEGG_LIMONENE_AND_PINENE_DEGRADATION- and KEGG_BUTANOATE_METABOLISM-associated pathways than the control group (Fig. 3G).
Enrichment analyses (Gene Ontology [GO]/Kyoto Encyclopedia of Genes and Genomes [KEGG])
To investigate the biological functions of the most significant module genes associated with WGCNA and glutamine metabolism compared with differential genes in atherosclerosis and normal controls, functional enrichment analysis was performed. GO results show, The most significant WGCNA and glutamine metabolism related module gene muscle contraction, muscle system process, regulation of actin filament—based process (BP), contractile fiber, myofibril, sarcomere(CC), actin binding, actin filament binding, integrin binding(MF) enrichment (Fig. 4A); The enriched KEGG pathways included Focal adhesion, Vascular smooth muscle contraction, Regulation of actin cytoskeleton, etc. (Fig. 4B). The GO results of differential genes in atherosclerosis and normal control showed that the genes were leukocyte cell–cell adhesion, T cell activation, leukocyte migration(BP), tertiary granule, secretory granule membrane, specific granule(CC), actin binding, immune receptor activity, integrin binding(MF) enrichment (Fig. 4C); Enriched KEGG pathways include Rheumatoid arthritis, Chemokine signaling pathway, Cell adhesion molecules and so on (Fig. 4D).
To investigate the biological functions of GLN-associated DEGs, GO term along with KEGG pathway enrichment analyses were performed. According to results of the GO term analysis, the genes exhibited strong enrichment in muscle contraction, muscle system process, and muscle tissue development (biological process [BP]); contractile fiber, myofibril, and sarcomere (cellular component [CC]); and actin binding, structural constituent of muscle, and transmembrane transporter binding (molecular function [MF]) (Fig. 5A–D).
The enriched KEGG pathways included dilated cardiomyopathy (DCM), hypertrophic cardiomyopathy (HCM), and cAMP signaling pathway (Fig. 5E).
Protein–protein interaction (PPI) network analysis and hub gene screening
To comprehend the interactions among the GLN-associated DEGs, a PPI network was built. Twelve approaches were adopted to elucidate highly correlated genes in the PPI network. Based on the intersections of the top 100 genes, obtained using all 12 methods, we identified 27 hub genes: ACTN2, TPM2, FLNC, MYH11, ITGA7, DMD, PAK3, LMOD1, GNAI1, MYLK, PLN, CTPS1, CNN1, MYH10, NLGN1, ADAMTSL3, MYOCD, MICU3, PPP1R14A, MAP2, OGN, PDE8B, RGS5, MAP1B, ITGA9, AKAP6, and SMTN (Fig. 6).
Hub gene validation
To further verify the diagnostic value of hub gene, ROC curve was used to verify hub gene, and FLNC (AUC = 0.8287), AKAP6 (AUC = 0.8139), LMOD1 (AUC = 0.8236), DMD (AUC = 0.8394) were found. ACTN2 (AUC = 0.8588), GNAI1 (AUC = 0.8222), CTPS1 (AUC = 0.8368), MAP1B (AUC = 0.8241), ITGA7 (AUC = 0.7965), ITGA7 (AUC = 0.7965), the area under ROC curve (AUC) values of ADAMTSL3 (AUC = 0.8218) and ITGA9 (AUC = 0.8014) were both greater than 0.6 (Fig. 7A). We used tenfold cross-examination to verify the diagnostic efficacy of the AUC model, and found that the AUC value was 0.957 (Fig. 7B), indicating that the AUC model has good diagnostic efficacy, which indicates that the hub gene has the differential ability as a potential biomarker of atherosclerosis. Among them, CNN1 (AUC = 0.874 (0.801, 0.947)) has good diagnostic value. Next, we used delong test to detect whether there were statistical differences in ROC curves between CNN1 and other hub genes, and the results showed that most hub genes were significantly different from CNN1(p < 0.05), suggesting that CNN1 might be a biomarkers of atherosclerosis (Fig. 7C,D, Supplementary Fig. S1A–M).
Immune cell infiltration (ICI)
Considering that ICI may play a crucial role in the pathogenesis of atherosclerosis, we examined the correlations between atherosclerosis/control samples and infiltrated immune cells. The atherosclerosis samples exhibited a markedly higher degree of infiltration than healthy controls, for most immune cells (Fig. 7A). A positive correlation was observed between the majority of immune cells (Fig. 8B). In addition, each hub gene showed a marked correlation with the corresponding immune cells (Fig. 8C–E). Notably, significant correlations were observed between MYLK and CD56bright natural killer (NK) cells (R = 0.881, p < 0.001), SMTN and CD56bright NK cells (R = − 0.889, p < 0.001), and MYLK and Th1 (R = − 0.898, p < 0.001) (Fig. 8C–E).
Signaling pathways participated in signature genes
We used GSVA to examine the differences in 50 HALLMARK signaling pathways in atherosclerosis and control samples.
In the atherosclerosis samples, the expression of 31 HALLMARK signaling pathways was markedly upregulated: HALLMARK_ADIPOGENESIS, HALLMARK_ALLOGRAFT_REJECTION, HALLMARK_ANGIOGENESIS, HALLMARK_APICAL_SURFACE, HALLMARK_APOPTOSIS, HALLMARK_COAGULATION, HALLMARK_COMPLEMENT, HALLMARK_DNA_REPAIR, HALLMARK_E2F_TARGETS, HALLMARK_ESTROGEN_RESPONSE_LATE, HALLMARK_G2M_CHECKPOINT, HALLMARK_GLYCOLYSIS, HALLMARK_HEME_METABOLISM, HALLMARK_HYPOXIA, HALLMARK_IL2_STAT5_SIGNALING, HALLMARK_IL6_JAK_STAT3_SIGNALING, HALLMARK_INFLAMMATORY_RESPONSE, HALLMARK_INTERFERON_ALPHA_RESPONSE, HALLMARK_INTERFERON_GAMMA_RESPONSE, HALLMARK_KRAS_SIGNALING_UP, HALLMARK_MTORC1_SIGNALING, HALLMARK_MYC_TARGETS_V2, HALLMARK_NOTCH_SIGNALING, HALLMARK_P53_PATHWAY, HALLMARK_PEROXISOME, HALLMARK_PI3K_AKT_MTOR_SIGNALING, HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY, HALLMARK_TNFA_SIGNALING_VIA_NFKB, HALLMARK_UNFOLDED_PROTEIN_RESPONSE, HALLMARK_UV_RESPONSE_UP, HALLMARK_XENOBIOTIC_METABOLISM. Conversely, the expression of the HALLMARK_MYOGENESIS pathway was markedly downregulated in the atherosclerosis samples (Fig. 9A).
We also analyzed the correlations between the top 5 most-significant differentially expressed hub genes and 50 HALLMARK signaling pathways. ACTN2 was associated with multiple pathways, including HALLMARK_UV_RESPONSE_UP and HALLMARK_UV_RESPONSE_UP (Fig. 9B).
Interaction analysis on hub genes
We created a PPI network for the signature genes using the GeneMANIA database and identified 27 genes (Fig. 10A). To further investigate the function of signature genes, GO and KEGG analyses were performed on 47 genes comprising 27 hub and 20 associated genes. According to the GO results, the genes were strongly enriched in regulation of axonogenesis, regulation of cell shape, and regulation of vascular smooth muscle (VSM) cell migration (BP); cell division site, cortical actin cytoskeleton, and caveola (CC); and cyclic-nucleotide phosphodiesterase activity, transcription coactivator activity, and transmembrane transporter binding (MF) (Fig. 10B). KEGG analysis revealed that DCM, VSM contraction, HCM, focal adhesion, regulation of actin cytoskeleton, arrhythmogenic right ventricular cardiomyopathy, adrenergic signaling in cardiomyocytes, cGMP-PKG signaling pathway, and oxytocin signaling pathway were the main enriched pathways (Fig. 10C).
Validation of the gene set
The expression levels of the key genes validated enrichment for gene set muscle tissue development, muscle system process, and muscle contraction. The atherosclerosis samples displayed low expression levels for most of the key genes (Fig. 11).