The clinical information and RNA-seq data of patients with IA were downloaded from the GEO (http://www.ncbi.nlm.nih.gov/geo) database [32]. After filtering, five GEO cohorts, including GSE75436 (15 IA and 15 control samples), GSE26969 (3 IA and 3 control samples), GSE13353 (19 IA samples), GSE15629 (14 IA and 5 control samples), and GSE54083 (13 IA and 10 control samples, Agilent platform), were included in this study.

In the GEO cohorts, only GSE54083 was collected from the Agilent platform, while the other four cohorts came from the Affymetrix platform. Thus, the expression matrix of the Affymetrix cohorts was subjected to background correction, quantile normalization, and log2 transformation to eliminate differences between groups using the robust multi-array average (RMA) algorithm in the “oligo” R package [33]. Meanwhile, the Agilent cohort was also performed with the quantile normalization and log2 transformation using the same method. Then, the probes were mapped to the genes based on the annotation files. The median value was taken to represent the gene expression level when one gene matched multiple probes, while the probe was removed when it mapped numerous genes. Further, the Combat algorithm was employed to eliminate the batch effects between different datasets [33]. Finally, 64 IA (GEO-IA) and 33 control samples were included in this study.

The 13 oxidative stress-related pathways (GO accession was shown in Supplementary Materials) and the “hallmark gene sets” of ROS and inflammatory response pathways from The Molecular Signatures Database (MsigDB, http://www.gsea-msigdb.org/gsea/downloads.jsp) were downloaded to assess the activation of oxidative stress pathway in IA samples [34]. In addition, the gene lists of amino acid metabolism pathways, carbohydrate metabolism, and lipid metabolism were downloaded from the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/) database using the keggGet function in the “KEGGREST” R package [35]. Then, we used the “GSVA” R package to perform the single-sample gene set enrichment analysis (ssGSEA) for calculating pathway enrichment score and for Gene Ontology (GO) and KEGG analyses [36].

The median absolute deviation (MAD) of encoding protein genes in each sample was calculated to obtain the top 70% of genes with potential expression dysregulation. These genes were subjected to WGCNA by using the “WGCNA” R package [37] to calculate the proper soft threshold (β) to ensure that the network was scale-free [38]. In addition, the topological overlap matrix (TOM) was converted from the weighted adjacency matrix to reduce noise and false correlation [39], and the Dynamic Tree Cut method was used to section the modules further [40]. Finally, the correlation between the modules and the ROS score was analyzed to determine the oxidative stress-related gene modules.

Under the thresholds of adj. p 0.05, false discovery rate (FDR) 0.05, and $|$log2FC$|$ $>$1, the differential expressed gene (DEG) analysis on the merged expression profile was performed in the “limma” R package [41]. Subsequently, the GO, KEGG, and GSEA enrichment analyses were conducted to analyze the functions of the DEGs [36, 42].

The random forest (RF) algorithm was run in the “randomForest” R package [43]. Lasso logistic regression was performed in the “glmnet” R package with minimal lambda as the optimal [44]. Additionally, the SVM classifier was developed by the “e1071” R package, and then SVM-RFE was performed using the RFE function in the “caret” R package [45] based on 10-fold cross-validation [46]. The timeROC R package [47] was applied to perform the receiver operating characteristic (ROC) analysis with the area under curve (AUC).

The ESTIMATE algorithm, including the stromal score, immune score, and ESTIMATE score, was used to assess the immune infiltration using the “estimate” R package [38]. The CIBERPSORT algorithm was used to evaluate the infiltration of 22 immune cells in the “CIBERPSORT” R package [48].

Human brain vascular smooth muscle cells (VSMCs) from the human brain were purchased from Bnbio (BNCC102172, Beijing, China). These cells were grown in Smooth Muscle Cell Medium (SMCM, Sciencell, San Diego, CA, USA) supplemented with 2% fetal bovine serum (FBS, Sciencell, San Diego, CA, USA), 5 mL of smooth muscle cell growth factor (SMCGS, Sciencell, San Diego, CA, USA), and 5 mL of penicillin/streptomycin (P/S, Sciencell, San Diego, CA, USA) in a 5% CO₂ incubator at 37 °C. Cells have been STR identified and the mycoplasma detection results are negative.

The TRIzol reagent extraction kit (Invitrogen, Thermo Fisher Scientific, Inc., Carlsbad, CA, USA) was utilized to isolate total RNA from both human brain VSMCs and control cells (untreated group). RNA concentrations were determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Carlsbad, CA, USA). The conditions for the PCR reaction were set as follows: an initial step at 95 °C for 30 seconds (s), followed by at 95 °C for 10 s, at 60 °C for 30 s, and at 70 °C for 34 s for a total of 40 cycles. All the data were calculated using 2^{-Δ⁢Δ⁢CT} and normalized to GAPDH. The experiments were performed in triplicate. The primers used are listed in Table 1.

The R software (v3.6.3, New York, NY, USA) package was used in all statistical analyses and data visualization in this study. Spearman’s method and the Wilcoxon test were adopted for correlation analysis and assessing significant differences (p 0.05 was defined as statistically significant). SangerBox (http://sangerbox.com/home.html) was also used to provide some technical support.

To explore the potential mechanisms of oxidative stress pathways in IA, we found that compared to the normal samples, the activation of oxidative stress pathways (Fig. 1A), such as the neuronal apoptosis signaling pathway, cellular regulation pathway, and regulation of response in response to oxidative stress was significantly increased in patients with IA. Similar results were also observed in the ssGSEA score (Fig. 1B), as the neuronal apoptosis signaling pathway and other pathways score was noticeably higher in IA samples than in control samples. The hallmark score of the ROS pathway in IA samples was also significantly higher than that in the control samples (Fig. 1C). Subsequently, we analyzed the correlation between the 13 oxidative stress pathways in IA samples and control samples and found that these pathways were closely correlated with each other in the two sample types (Fig. 1D). Specifically, the correlation coefficient between the regulation of cellular response to oxidative stress pathway and regulation of response to oxidative stress pathway in IA samples (lower triangle) and control samples (upper triangle) was 0.96 and 0.94, respectively. These results indicated that the abnormal activation of oxidative stress pathways may be closely associated with the development of IA.

Fig. 1.

The enrichment analysis of oxidative stress pathways in intracranial aneurysm (IA) progression. (A) The single-sample gene set enrichment analysis (ssGSEA) of 13 oxidative stress pathways in Gene Expression Omnibus (GEO)–IA cohort. (B) The difference in ssGSEA between IA and control samples in 13 oxidative stress pathways. (C) The enrichment of hallmark reactive oxygen species (ROS) pathway in IA and control samples. (D) Correlation analysis between 13 oxidative stress pathways in IA and control samples (upper triangle represented the control samples; lower triangle represented the IA samples). ns p $>$ 0.05, **p 0.01, ***p 0.001, ****p 0.0001.

The ssGSEA score of metabolism pathways was computed to explore the metabolic difference. In carbohydrate metabolism, multiple metabolic pathways, such as the fructose and mannose and amino sugar and nucleotide sugar metabolism, starch and sucrose, galactose, pentose phosphate pathway (PPP), were significantly activated in the IA samples, while the pyruvate, propanoate, ascorbate and aldarate metabolism pathways butanoate, glyoxylate and dicarboxylate were significantly activated in the control samples (p 0.05, Fig. 2A). In lipid metabolism, the sphingolipid, glycerophospholipid and ether lipid metabolism pathway were significantly activated in the IA samples, while the linoleic acid, alpha-linolenic acid metabolism and fatty acid degradation pathway were significantly activated in the control samples (p 0.05, Fig. 2A). In amino acid metabolism, most pathways were significantly inhibited in the IA samples, such as glycine, serine, threonine, histidine, phenylalanine metabolism, only leucine and isoleucine biosynthesis, the cysteine and methionine metabolism and valine were significantly activated (p 0.05, Fig. 2A). After visualizing the scores of these metabolism pathways via a heatmap, most amino acid metabolism pathways were significantly inhibited in the IA samples, while the carbohydrate metabolism pattern was significantly different in the IA samples from that in control samples (p 0.05, Fig. 2B). Additionally, we investigated the correlation between the metabolism and oxidative stress pathways, and the results showed that most carbohydrate and lipid metabolism pathways were significantly positively correlated with the oxidative stress pathways, whereas most amino acid metabolism pathways were significantly negatively correlated with the oxidative stress pathways (p 0.05, Fig. 2C).

Fig. 2.

The enrichment analysis of metabolism pathways in IA progression. (A) The enrichment of carbohydrate, lipid, and amino acid metabolic pathways in IA and control samples. (B) The carbohydrate, lipid, and amino acid metabolic enrichment heatmap for IA and control samples. (C) The correlation analysis between oxidative stress pathways and metabolism pathways in IA and control samples (upper triangle represented the control samples; lower triangle represented the IA samples). ROS, reactive oxygen species; OSR, oxidative stress response. ns p $>$ 0.05, *p 0.05, **p 0.01, ***p 0.001, ****p 0.0001.

Furthermore, we calculated the immune infiltration-related score to explore differences in tumor microenvironment (TME) difference in the two types of samples. The ESTIMATE analysis showed that IA samples had significantly higher stromal, immune, and ESTIMATE scores than the control samples (Fig. 3A); meanwhile, the inflammation response pathways, including the MAPK, PI3K–Akt–mTOR, and IL6–JAK–STAT3 signaling pathways, also presented a significantly high ssGSEA score in the IA samples (Fig. 3B). In addition, CD8 T cells, gamma delta T cells, macrophages M2 and activated mast cells showed noticeably high infiltration levels, contributing to an immune response in the IA samples, while the naïve B cells, resting mast cells, and resting CD4 memory T cells had been significantly infiltrated in the control samples (Fig. 3C). The correlation analysis between the immune infiltration and oxidative stress pathways demonstrated that these inflammation-related pathways were closely correlated with theses oxidative stress pathways (Fig. 3D).

Fig. 3.

Immune infiltration difference between IA and control samples. (A) ESTIMATE analysis between IA and control samples. (B) The ssGSEA of the inflammation response pathway in IA and control samples. (C) CIBERPSORT analysis between IA and control samples. (D) The correlation analysis between oxidative stress pathways and immune infiltration in IA and control samples (upper triangle represented the control samples; lower triangle represented the IA samples). ns p $>$ 0.05, *p 0.05, ** p 0.01, ***p 0.001, ****p 0.0001.

We performed WGCNA to classify oxidative stress-related gene modules. The top 70% of genes above the MAD were subjected to WGCNA under the soft threshold $\beta$ of 6 (Supplementary Fig. 1A,B) to ensure the scale-free nature of the network. Finally, 17 gene modules were sectioned (Supplementary Fig. 1C,D). The correlation between the gene modules and other clinical features was calculated (Fig. 4A). As the blue module exhibited a significantly positive correlation with the sample type and ROS hallmark pathway (p 0.001, r = 0.7), it was considered as the module being closely associated with oxidative stress. In addition, the correlation coefficient between module membership (MM) and gene significance (GS) in the blue module reached 0.83 (p 0.001, Fig. 4B), suggesting that these genes were closely related to the ROS phenotype and the gene module had a high quality. Further, we screened 610 key genes from the blue module (MM $>$0.7 and GS $>$0.4) and performed GO and KEGG enrichment analyses. The KEGG analysis demonstrated that these genes were closely associated with the immune pathways of the lysosome, phagosome, T cell killing, and chemokine signaling activation (Fig. 4C). These genes were related to the immunity and oxidative stress pathway activation of neutrophils and T cells in the biological process in the GO analysis (Fig. 4D). The cell component in the GO indicated that these genes were associated with the vesicle lumen, membrane, granule, lysosome lumen pathways (Fig. 4E). These genes were close with receptors, cytokine binding, and enzymatic activity in the molecular function of the GO analysis (Fig. 4F). Thus, these oxidative stress-related genes were used for further study.

Fig. 4.

WGCNA for oxidative stress-related gene module. (A) Correlation analysis between module eigengenes and clinical traits in the Gene Expression Omnibus (GEO). (B) The high correlation between GS and MM in the blue module (p 0.0001) in the GEO. (C) Kyoto encyclopedia of genes and genomes (KEGG) analysis of oxidative stress-related genes. (D) Biological process of GO analysis of oxidative stress-related genes. (E) Cell component in the GO analysis of oxidative stress-related genes. (F) Molecular function in the GO analysis of oxidative stress-related genes.

Analysis of the DEGs between the IA and control samples determined 770 DEGs (Fig. 5A), among which 380 genes were upregulated and 390 were downregulated (Fig. 5B). The function enrichment of these upregulated genes indicated that they were closely associated with the neutrophil immune activation, cytokine production, biotic stimulus, and oxidative stress pathways (Fig. 5C). Interestingly, neutrophil peroxidases have been reported as a potent source of oxidative stress [49]. These downregulated genes were closely associated with muscle and circulation regulation and the tissue differentiation pathway (Fig. 5D). The results of the GSEA showed that the oxidative stress pathways, such as the neuron intrinsic apoptotic signaling pathway in the intrinsic apoptotic signaling pathway, cellular response to oxidative stress, response oxidative stress in response to oxidative stress, were significantly enriched in the IA samples (Fig. 5E). Overall, these results suggest that the abnormal activation of the oxidative stress pathway was significantly associated with IA progression.

Fig. 5.

Differential expression analysis in IA and control samples. (A) Identification of differentially expressed genes (DEGs) between IA and control samples. (B) The heatmap of DEGs in IA and control samples. (C) The enrichment analysis of upregulated genes. (D) The enrichment analysis of downregulated genes. (E) The GSEA of IA samples.

We obtained 209 upregulated key oxidative stress-related genes from the intersection between the 171 upregulated DEGs and the 610 key oxidative stress-related genes identified by WGCNA (Fig. 6A). Thus, we aimed to identify the key biomarkers of oxidative stress for IA from these 209 genes. First, 13 potential biomarker candidate genes were detected by running the SVM-RFE algorithm (Fig. 6B). Next, we used the Lasso–Cox regression analysis to identify 10 potential biomarker candidate genes (Fig. 6C). Finally, the RF algorithm screened 10 crucial genes associated with the oxidative stress response (Fig. 6D). The Venn plot of these above candidate genes showed that FLVCR2, SDSL, TBC1D2, and SLC31A1 in the intersection (Fig. 7A) could be defined as the key markers associated with oxidative stress for the diagnosis of IA.

Fig. 6.

Screening key oxidative stress-related biomarkers for IA patients. (A) Venn plot of the upregulated DEGs and key oxidative stress-related module genes. (B) SVM-RFE-screened biomarkers. (C) Lasso-logistic regression algorithm for screening diagnostic markers. (D) Using random forest (RF) algorithm for screening biomarkers.

Fig. 7.

Verification of the diagnostic biomarkers. (A) Venn plot of key biomarkers identified by SVM-RFE, Lasso, and RF algorithms. (B) ROC analysis of the four biomarkers for verifying diagnostic performance. (C) Receiver operating characteristic (ROC) analysis of the combination of the four biomarkers for verifying diagnostic performance. (D) The expression differences in the four biomarkers in the IA and control samples. ****p 0.0001.

We performed the ROC analysis for each gene to explore the clinical value of these four genes in diagnosing IA. The results showed that the AUC values of FLVCR2, SDSL, TBC1D2, and SLC31A1 were as high as 0.901, 0.913, 0.921, and 0.911, respectively (Fig. 7B), indicating that these four biomarkers had a high predictive value accuracy. Moreover, the combination of these four genes reached an AUC of 0.951 (Fig. 7C), which showed that the combined use of these four biomarkers also had excellent classification performance and clinical diagnostic value. In addition, these genes were significantly highly expressed in IA patients (Fig. 7D). Vascular smooth muscle cells play an important role in the formation and development of IA [50]. Therefore, to validate the potential impact of these four biomarkers on IA development further, we detected their expression in human brain VSMCs using qRT-PCR. As shown in Supplementary Fig. 2, the expressions of FLVCR2, SDSL, TBC1D2, and SLC31A1 were significantly upregulated in human brain VSMCs compared to the controls.

The correlation analysis between the four biomarkers and oxidative stress pathways was performed. The results showed that the FLVCR2 expression was significantly positively correlated with multiple oxidative stress pathways, such as the regulation of the response, cellular response, and intrinsic apoptotic signaling pathway in response to oxidative stress (p 0.0001, Fig. 8A). Meanwhile, the SDSL expression was significantly positively correlated with the neuron intrinsic apoptotic signaling pathways and the cellular and regulatory responses to oxidative stress (p 0.0001, Fig. 8B). The expression levels of TBC1D2 (Fig. 8C) and SLC31A1 (Fig. 8D) were closely related to multiple oxidative stress pathways. In addition, the correlation analysis between the ROS pathway and these biomarkers showed a significant positive correlation relationship, with the correlation coefficients for FLVCR2, SDSL, TBC1D2, and SLC31A1 of 0.609, 0.514, 0.44, and 0.445, respectively (Fig. 8E).

Fig. 8.

The relationship between the four biomarkers and oxidative stress activation. (A) The correlation between oxidative stress pathways and FLVCR2 expression. (B) The correlation between oxidative stress pathways and SDSL expression. (C) The correlation between oxidative stress pathways and TBC1D2 expression. (D) The correlation between oxidative stress pathways and SLC31A1 expression. (E) The correlation between ssGSEA of ROS pathway and the expression of the four biomarkers.

Subsequently, the relationship between the biomarkers and immune infiltration [51] was analyzed, and we observed that the expression of FLVCR2 was significantly positively correlated with the infiltration of activated CD4 memory T cells, macrophage M0 and M2 but significantly negatively correlated with activated NK cells and naïve B cells (p 0.01, Fig. 9A). Meanwhile, the expression of SDSL was significantly positively correlated with the infiltration of macrophage M0 and M2 but significantly negatively correlated with activated NK cells (p 0.01, Fig. 9B). Similar results also were observed for TBC1D2 (Fig. 9C) and SLC31A1 (Fig. 9D), as we found that the expression of these two genes was significantly positively correlated with the infiltration of macrophage M0 and M2 but significantly negatively correlated with activated NK cells. These immune cells contributed to a suppressive TME of IA to promote cancer progression.

Fig. 9.

The relationship between the four biomarkers and immune infiltration. (A) The correlation between immune cell infiltration and FLVCR2 expression. (B) The correlation between immune cell infiltration and SDSL expression. (C) The correlation between immune cell infiltration and TBC1D2 expression. (D) The correlation between immune cell infiltration and SLC31A1 expression.

Correlation analysis between the biomarkers and metabolism pathways revealed that these four biomarkers were closely associated with multiple metabolism pathways (Fig. 10). For example, the expression of FLVCR2 and SDSL was significantly positively correlated with sphingolipid, glycerophospholipid, glycerolipid, fructose, and mannose metabolism, amino sugars and nucleotide sugars, ether lipids, arachidonic acid, linoleic acid, and galactose but they were all significantly negatively correlated with several amino acid metabolic pathways, such as glycine, serine, threonine, and tyrosine metabolism, and the valine, leucine, and isoleucine degradation pathways (Fig. 10). The expression of TBC1D2 and SLC31A1 was significantly positively correlated with the galactose metabolic pathway, linoleic acid, and amino sugars and nucleotide sugars (Fig. 10). Interestingly, the metabolism of linoleic acid, glycolysis pathways, fructose and mannose metabolism, amino sugars and nucleotide sugars, galactose, and pentose phosphate (PPP) were closely associated with the expression of these four biomarkers, indicating that these pathways played critical roles in promoting IA progression.

Fig. 10.

The relationship between the four biomarkers and metabolic pathways.