The overall workflow of this study is depicted in Figure 1 . Our study retrospectively recruited seven independent IA cohorts from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database, including GSE122897 (21 ruptured intracranial aneurysm (RIA) tissues, 21 UIA tissues, and 16 control cerebral artery tissues), GSE54083 (13 IA tissues, and 10 control cerebral artery tissues), GSE75436 (15 IA tissues, and 15 control cerebral artery tissues), GSE13353 (11 RIA tissues and 8 UIA tissues), GSE66239 (7 IA tissues and 10 control cerebral artery tissues), GSE36791 (43 blood samples from patients with RIA, and 18 blood samples from control subjects without IA), and GSE159610 (25 blood samples from patients with UIA, and 22 blood samples from control subjects without IA). The detailed baseline information was summarized in Supplementary Table S1 . Among these cohorts, the GSE122897 dataset was utilized for deciphering immune-relevant alterations in IA local vessels, and GSE54083, GSE75436, and GSE13353 datasets were employed for verification. The precious miRNA array dataset GSE66239 was applied for identifying deregulated microRNA (miRNAs) in IA. A total of 108 blood samples generated by integrating the GSE36791 and GSE159610 datasets were denoted as meta-cohort and utilized to develop and validate a blood-based PFDGS model for IA risk assessment.

For GSE122897, RNA read count data were normalized by variance stable transformation (VST) and further differential expression analyzed using the DESeq2 R package (24). Genes with an absolute value of log2FoldChange >0.5 and false discovery rate (FDR) <0.05 were considered as differentially expressed genes (DEGs). The microarray data of GSE54083, GSE75436, GSE13353, GSE66239, and GSE36791 datasets were preprocessed via quantile normalization and log2 transformation. The RNA-seq raw read counts in GSE159610 were converted to transcripts per kilobase million (TPM) values and further log-2 transformed, which is more similar to those resulting from microarrays and more comparable between samples. The ComBat algorithm implemented in the sva R package (25) was applied to remove batch effects from non-biological technical biases between GSE36791 and GSE159610 datasets. The principal component analysis (PCA) before and after batch correction was shown in Supplementary Figure S1 . Ultimately, a total of 14,070 genes were retained in the meta-cohort for subsequent analysis in this study. The GENCODE database (https://www.gencodegenes.org/) was utilized for mRNA and long non-coding RNA (lncRNA) annotations. Gene differential expression analysis of GSE66239 and meta-cohort were performed using limma R package (26).

According to the expression level of immune cell-specific marker genes (27), the relative abundance of 28 immune cells in each sample was quantified by single-sample gene set enrichment analysis (ssGSEA) algorithm implemented in the GSVA R package (28), which is broadly utilized in immune infiltration-relevant bioinformatics studies (13–15, 29). Please refer to the Supplementary Table S2 for details of the gene sets marking 28 immune cells.

To explore the expression patterns of immune-related molecules in IA, we searched the ImmPort database (https://www.immport.org/) and obtained 1775 immune-related genes ( Supplementary Table 3 ) with regard to BCR signaling pathway, TCR signaling pathway, NK cell cytotoxicity, chemokines, cytokines, IFNs, ILs, TGF-b family members, TNF family members, corresponding receptors, and other immune response molecules. Immune-related genes overlapping with DEGs were defined as immune-related DEGs in IAs.

In order to ascertain the status of immune-related pathways in IAs, we performed gene set enrichment analysis (GSEA) based on 27 Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway gene sets, which were retrieved from the MSigDB (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). Detailed pathways and gene lists were depicted in Supplementary Table S4 . Immune-related DEGs and co-expression module genes were functionally annotated by Gene Ontology (GO) and KEGG enrichment analysis to understand their biological functions. The above procedures were performed using clusterProfiler and fgsea R packages (30). Permutations were set to 10,000 to obtain normalized enrichment scores (NES) in GSEA. Gene sets with FDR <0.05 were considered significantly enriched. Moreover, the GSEA algorithm was also employed to decipher potential biometric differences between high- and low-risk groups distinguished by PFDGS scores.

The STRING database (https://www.string-db.org/) was applied to generate a PPI network of the immune-related DEGs, which reflects the interaction relationship between the proteins encoded by these genes. Only those interaction pairs with a combined score >0.7 were picked up as remarkable. Hub genes in the PPI network were further identified utilizing the degree algorithm of the cytoHubba plug-in of Cytoscape software.

This study was approved by the Ethical Committee of Zhengzhou University People’s Hospital, China. Five IA tissue specimens and four superficial temporal artery (STA) controls were obtained from patients undergoing IA clipping resection in the Department of Neurosurgery, Zhengzhou University People’s Hospital. All patients were aged >18 years and gave written informed consent.

Paraffin-embedded specimen sections were incubated at 65 °C for 20 min and deparaffinized with xylene and alcohol before staining experiments. Hematoxylin and eosin (HE) stained sections were exploited to ascertain the gross structure of the tissues. For IHC, sections were sequentially incubated with 3% H₂O₂, goat serum, primary antibody, biotinylated second antibody, and horseradish enzyme-labeled streptavidin according to the IHC SP kit (ZSGB-BIO: OriGene Technologies, Cat. No. SP-9000) instructions, and then stained by diaminobenzidine (DAB) and hematoxylin. Primary antibodies used in the study were as follows: anti-CXCL10 (1:500; Proteintech, Cat. No. 10937-1-AP), anti-IL6 (1:500; Abcam, Cat. No. ab6672), anti-IL10 (1:500; Proteintech, Cat. No. 20850-1-AP), anti-STAT1 (1:200; Abcam, Cat. No. ab109320), anti-VEGFA (1:100; Abcam, Cat. No. ab52917). IHC results were quantified using the integrated optical density (IOD) method via Image-Pro Plus 6.0 software.

The upstream ceRNA regulatory networks of hub genes were constructed through the following pipelines: Firstly, miRNAs targeting five hub genes were predicted via miRWalk (http://mirwalk.umm.uni-heidelberg.de/) and filtered with miRWalk Score >0.95. Secondly, differently expressed miRNAs (DEmiRNAs) in IA were identified based on the miRNA array dataset (GSE66239) via limma R package with the criteria of abs(log2FoldChange) >0.5 and P-value <0.05. DEmiRNAs overlapping with the predicted miRNAs were considered regulatory miRNAs of hub genes and included in the subsequent analysis. Thirdly, the target lncRNAs of these miRNAs were further predicted according to the lncBase v.2. LncRNAs with abs(log2FoldChange) <0.5 and predicted binding pairs with miTG-score <0.95 were eliminated. Fourthly, the miRNA-mRNA and miRNA-lncRNA binding pairs obtained from the above steps were merged into multiple lncRNA-miRNA-mRNA regulatory axes. Only regulatory axes with significant positive correlations (R >0.3 and P <0.05) between lncRNAs and mRNAs were retained, and then, the ceRNA regulatory networks of the five hub genes were completely constructed.

The potential drugs were predicted via DGIdb (https://dgidb.genome.wustl.edu/), an open-source database that converged the known or potential interaction between drugs and genes. Only drug-gene pairs with established interaction types and supported by literature was selected and visualized by Cytoscape software.

The advanced WGCNA method was utilized to identify gene co-expression modules associated with immune cell infiltration. Based on the DEGs between IA and controls, gene co-expression networks were generated using WGCNA R package (31). Briefly, genes with similar expression patterns were assigned to a co-expression module based on a weighted correlation adjacency matrix and cluster analysis. An appropriate soft threshold β was calculated to meet the criteria for a scale-free network and construct a weighted adjacency matrix. Further, the weighted adjacency matrix was converted into a topological overlap matrix (TOM), and the corresponding dissimilarity was generated (1-TOM). The dynamic tree cutting approach was utilized for module recognition, and modules with less than 0.3 dissimilarities were merged. The relationship between module eigengenes values and immune cell abundance was evaluated via Pearson correlation. Modules that were significantly associated with the vast majority of immune cells were considered crucial immune-related modules and were selected for in-depth investigation.

Peripheral blood DEGs between IAs and controls were identified in the GSE36791 and GSE159610 combined meta-cohort, with limma-FDR less than 0.05 as the criteria. Then, peripheral blood DEGs were taken to intersect with dysregulated pathological feature-related genes, comprising upregulated immune-related genes and functional module genes recognized by WGCNA. Genes that presented consistent expression trends in peripheral blood and vascular tissue were incorporated into the model development. Moreover, we transformed gene expression into z-score in the meta-cohort prior to constructing the PFDGS model, which enhanced the comparability between different datasets.

The meta-cohort was randomly divided into a training set and validation set in a 7:3 ratio. The machine learning-based development procedure for the PFDGS model in the training set was as follows: (a) Logistic regression analysis was performed on the pathological feature-related genes, and the genes dramatically associated with IA rupture were selected for subsequent modeling. (b) The random forest algorithm was employed to further screen critical genes, bounded by relative importance greater than 0.5 for variable filtering. (c) Ultimately, the PFDGS signature was constructed using the least absolute shrinkage and selection operator (LASSO) regression based on critical genes recognized from the random forest filter. The PFDGS risk score for each patient was calculated as follows:

where i is the key genes obtained from the LASSO regression, Expr(i) is the relative expression value (z-score) of i, and Coef(i) is the regression coefficient of i.

The diagnostic and predictive value of critical genes and the PFDGS model were estimated in the training set, validation set, and meta-cohort, respectively.

All data processing, statistical analysis, and plotting were conducted in R 4.1.0 software. Significance was assessed via Student’s t-test or Wilcoxon rank-sum test for comparisons of two groups and Kruskal-Wallis test for comparisons of three or more groups. Correlations between two continuous variables were determined using Pearson correlations. Random forest and LASSO regression were implemented via randomForest and glmnet R packages (32), respectively. Receiver operating characteristic (ROC) curve analysis applied to evaluate the diagnostic or predictive power of variables was performed via pROC R package (33). All statistical tests were two-tailed and P <0.05 was considered as statistically significant.

Using the ssGSEA algorithm, we first measured the relative infiltration abundance of 28 immune cell subpopulations in each sample of the GSE122897 cohort. Overall, most immune cells presented high infiltration levels in IAs, particularly in RIAs ( Figures 2A, B and Supplementary Figures S2A, B ), which can be figuratively described as an “immuno-thermal” microenvironment. HE staining revealed the histological features of IAs, mainly comprising internal hemorrhage, structural disorganization, as well as pronounced accumulation of immune cells, which evidenced our inference ( Figure 2D ).

To further ascertain immune cells with stable differential infiltration levels, we additionally assessed immune infiltration in the GSE54083, GSE75436, and GSE13353 datasets ( Supplementary Figures S2C-E ). Compared with normal controls, IAs presented remarkable higher infiltration of 11 effector immune cells, including macrophage, activated dendritic cell (DC), natural killer (NK) cell, NK T cell, CD56+ NK cell, myeloid-derived suppressor cell (MDSC), activated CD4 T cell, activated CD8 T cell, gamma delta (γδ) T cell, regulatory T (Treg) cell, and Type 1 T helper (Th1) cell ( Supplementary Figure S2A, C and E , Supplementary Table S5 ). Moreover, two important inflammatory cells, mast cell, and neutrophil, were also relatively more abundant in IAs, although not stably significant in validation datasets. Further comparison of RIA and UIA detected that the infiltration levels of activated CD4 T cell activated DC, MDSC, γδ T cell, Treg cell, monocyte, and macrophage were elevated in RIAs relative to UIAs ( Supplementary Figures S2B, D and Supplementary Table S5 ). In addition, no significant differences were observed when contrasting between men and women, suggesting that gender was not an interfering factor ( Supplementary Figure S2F ).

The correlation between immune cells was visualized in the heatmap and chord diagram ( Supplementary Figure S3 and Figure 2C ). The vast majority of immune cells exhibited a strong positive correlation with each other. Strikingly, among the above 14 vital dysregulated immune cells, two immunosuppressive cells were synergistically elevated with eight innate immune cells and four adaptive immune cells ( Figure 2C ), suggesting a specific microenvironment symptom in IAs: co-existence of immune activation and immunosuppression.

In GSE122897, GSEA performed on 27 KEGG signaling pathways showed that immune relevant pathways, such as JAK-STAT signaling pathway, Toll-like receptor signaling pathway, T cell receptor signaling pathway, NOD-like receptor signaling pathway, and chemokine signaling pathway, were remarkably activated in IAs ( Figure 2G ); while smooth muscle contraction relevant pathways encompassing vascular smooth muscle contraction and calcium signaling pathway were pronounced inhibited in IAs. These aberrant signaling pathways characterized the pathological features of immune activation and contractile smooth muscle cell (SMC) dysfunction or loss in IAs and may serve as pathway-based targets candidate for molecular therapy.

A total of 3310 DEGs between IAs and normal controls were identified in GSE122897, including 2746 differentially expressed mRNAs (DEmRNAs, Figure 2E ), 552 DElncRNAs ( Figure 2F ), and 12 unannotated RNAs. Subsequently, intersecting these DEGs with 1775 immune-related genes collected from the ImmPort, we detected 146 upregulated and 99 downregulated immune-related DEGs in IAs ( Figure 3A ). GO and KEGG enrichment analysis indicated that biological processes with respect to proliferation and activation of various leukocytes, cell chemotaxis and migration, and regulation of cytokine production, as well as immune relevant pathways such as cytokine-cytokine receptor interaction and JAK-STAT signaling pathways, were remarkably enriched in upregulated immune-related DEGs ( Supplementary Figures S4A, C ); whereas the downregulated genes were dramatically associated with the regulation of growth and development ( Supplementary Figures S4B, D ). We thus speculated that those upregulated immune-related DEGs might be major effectors mediating immune activation in IAs and deserved further attention.

A PPI network reflecting the association of proteins encoded by these upregulated immune-related DEGs was generated via the STRING database ( Figure 3B ). Using the cytoHubba plug-in of Cytoscape software, we extracted the Top5 genes ranked by degree method from the PPI network as immune-related hub genes, comprising IL6, IL10, STAT1, CXCL10, and VEGFA. As well-known, the expression products of these five hub genes are essential immunological molecules involved in immune response regulation. Moreover, correlation analysis indicated that all five hub genes were positively correlated with the 14 vital dysregulated immune cells we highlighted above, implying the immune cell recruitment potential of these hub genes ( Figure 3C ). Notably, anti-inflammatory cytokine IL10 was pronounced upregulated in IA and strongly correlated with immune infiltration, which suggested a microenvironment feature similar to our cell-based inference that immunosuppression and immune activation coexisted in IAs. To further verify the differential expression of hub genes at the protein level, we performed IHC on paraffin sections from five IA and four STA control samples. Likewise, the protein expression level of CXCL10 (P <0.001), IL6 (P <0.001), IL10 (P =0.008), STAT1 (P =0.002), and VEGFA (P <0.001) were remarkably higher in IAs ( Figures 3D, E ).

In the GSE66239 dataset, a total of 471 DEmiRNAs were identified, including 194 upregulated and 277 downregulated miRNAs ( Supplementary Figure S5A ). Considering that all five hub genes were upregulated, we only extracted downregulated DEmiRNAs for subsequent analysis. The flowchart of ceRNA network construction is displayed in Supplementary Figure S5B . According to the integrated lncRNA-miRNA-mRNA analysis pipeline we proposed above, five ceRNA networks were constructed to decipher the potential post-transcriptional regulation mechanism of hub genes ( Figure 4A ). In each network, lncRNAs may promote the expression of hub genes by inhibiting the corresponding miRNAs. For instance, CDKN2B-AS1 may inhibit has-miR-6827-5p and thus upregulate the expression of IL6. Notably, LINC00960 might be able to affect both CXCL10 and STAT1 through multiple miRNAs, including a consensus has-miR6801-5p. Likewise, IL10 and VEGFA shared four lncRNA regulators encompassing LINC1001, CTD-2006K23.1, RP4-647J21.1, and CSNK1G2-AS1, which may be vital RNA agents or targets candidate for anti-inflammatory and proendothelialization molecular therapy. In addition, the latent drugs interacting with hub genes were also explored in the DGIdb database. The IL6 inhibitors Olokizumab and Siltuximab, and the VEGFA inhibitors Aflibercept, Bevacizumab, Ranibizumab, and Pegaptanib Sodium were predicted to serve as potential therapeutic agents for IAs ( Figure 4B ).

The 3310 DEGs between IA and control were included in the WGCNA analysis, and the soft power of β =18 (scale-free R^{2 =} 0.89) was determined as soft-thresholding to acquire DEG co-expressed gene modules ( Figure 5A ). Ultimately, 739,394,179,509 and 268 DEGs were clustered into turquoise, brown, green, blue, and yellow modules, respectively; while the other 1221 genes were assigned to the grey region due to failed clustering. ( Figure 5B ). Among these gene modules, yellow module was significantly positively associated with aneurysm rupture (r =0.3, P =0.02) and nine vital immune cells, such as γδ T cell (r =0.37, P =0.004) and macrophage (r =0.32, P =0.01, Figure 5C ). Brown module was remarkably negatively correlated with ten vital immune cells such as CD56+ NK cell (r = -0.37, P =0.004) and macrophage (r = -0.3, P =0.02), and this module also presented a potential negative correlation with rupture status, although the P-value is not perfectly significant (r = -0.24, P =0.06). Taken together, these two modules were positively or negatively relevant to the most immune cells and were considered important transcriptionally dysregulated gene modules with regard to the altered immune microenvironment.

To further ascertain the biological function of the modules, GO and KEGG enrichment analyses were conducted for genes in the brown and yellow modules, respectively. Brown module genes were pronounced enriched in muscle contraction, regulation of cytosolic calcium ion concentration, calcium ion transmembrane import into the cytosol, cAMP signaling pathway, calcium signaling pathway, and so on, which can presumably be summarized as SMC contraction-related gene module ( Figures 5D, E ). Likewise, the yellow module was characterized as ECM metabolism-related gene modules, which were mainly enriched in the extracellular matrix organization, external encapsulating structure organization, collagen metabolic process, PI3K-Akt signaling pathway, ECM-receptor interaction, etc. ( Figures 5F, G ). Considering that yellow and brown modules contained mainly up- and down-regulated DEGs, respectively, we speculated that the altered immune microenvironment may accompany and even mediate two other essential pathologies in IA: contractile SMC dysfunction or loss, and ECM degradation and remodeling ( Figure 5H ).

Taken together, upregulated immune-related DEGs, contractile SMC loss-related DEGs, and ECM metabolism-related DEGs were considered essential pathological feature-related genes in IA, which may serve as promising biomarkers. To explore the diagnostic and predictive value of these genes and enhance their clinical applicability, we attempted to pursue shared DEGs in peripheral blood as candidate biomarkers. As a result, 27 DEGs with consistent alteration trends in blood and tissue were detected in the meta-cohort. Subsequently, the expression profiles of these 27 DEGs were subjected to the machine learning-based procedure to develop a pathological feature-derived gene signature (PFDGS). Univariate Logistic regression analysis of the 27 DEGs revealed 15 marker genes that were also significantly associated with IA rupture ( Figure 6A ), followed by random forest, screening out nine genes with relative importance greater than 0.5 ( Figures 6B, C ). The nine genes were incorporated into LASSO regression to generate the final PFDGS model. Based on the 5-fold cross-validation, the model reached an optimum when lambda was equal to 0.037, containing five key gene variables ( Figures 6D, E ). A risk score for each patient was calculated based on the expression of the five model genes with detailed formula as follow: PFDGS = 0.301*IL1R2 + 0.018*NFKBIZ + 0.277*S100A12 + 0.410*STK17B + 0.495*TPST2 + 0.772 ( Figure 6F ).

Patients with IA presented higher PFDGS scores than IA-free controls ( Figure 6G ); and when compared between RIAs and UIAs, PFDGS was relatively higher in RIAs ( Figure 6H ). To quantitatively assess the diagnostic and predictive value of the model, we conducted an ROC curve analysis of PFDGS and 15 initial candidate marker genes. When IA patients were set as positive and the controls as negative, the AUC of PFDGS in the training cohort (n =76), validation cohort (n =32), and meta-cohort (n =108) were 0.810, 0.838, and 0.813, respectively, indicating that the PFDGS had a good performance in the diagnosis of IA ( Figures 7A-C ). Likewise, when patients with RIA were set as positive and patients with UIA as negative, the AUC of PFDGS in the training cohort (n =48), validation cohort (n =20), and meta-cohort (n =68) were 0.909, 0.920, and 0.913, respectively, suggesting that the PFDGS also possessed a great performance in predicting IA rupture ( Figure 7D ). Moreover, PFDGS also displayed superior AUC values relative to 15 initial candidate marker genes alone for diagnosis or prediction, which implied the superiority of the gene-combined signature ( Figures 7E, F ). Overall, patients with high PFDGS scores may possess a higher risk of IA onset or rupture, requiring more vigilant attention or timely interventional treatment.

To further investigate the potential biological mechanisms underlying the PFDGS, we divided all samples in the meta-cohort into high- and low-risk groups according to the median PFDGS score and performed GSEA analysis. The high-risk group enriched numerous pathways associated with immune-inflammatory activation and metabolic disorders, such as leukocyte activation involved in immune response, macrophage activation, neutrophil extracellular trap formation, positive regulation of interleukin-1 production, NOD-like receptor signaling pathway, TNF signaling pathway, reactive oxygen species metabolic process, and regulation of actin cytoskeleton ( Figures 8A, B ). We subsequently measured immune cell abundance based on blood gene expression profile using ssGSEA, and detected 19 differentially infiltrated immune cells between high- and low- risk groups ( Figure 8C ). Of note, the high-risk group presented more abundant myeloid cells (e.g., neutrophils, eosinophils, and macrophages) but few lymphocytes (e.g., activated CD8 T cells, immature B cell, and NK T cell) relative to the low-risk group. This is similar to the special phenomenon of increased inflammatory cells but decreased lymphocytes in the peripheral blood of IA patients discovered in recent years. PFDGS was positively or negatively correlated with up- and down-regulated immune cells, respectively ( Figure 8D ), implying an ability of PFDGS to reflect the special alteration of immune cells in the peripheral blood of IA patients. In addition, we also investigated the association of PFDGS with matrix metalloproteinases (MMPs) and tissue inhibitors of metalloproteinases (TIMPs), two essential molecular clusters involved in ECM metabolism (34). In the meta-cohort, five genes encoding MMPs and TIMPs were detected, and MMP9 and TIMP2 were significantly upregulated in the high-risk group ( Figure 8E ). Correlation analysis indicated that PFDGS was positively associated with the expression of MMP9 (r =0.77, P <0.001) and TIMP2 (r =0.60, P <0.001; Figures 8F, G ), suggesting the potential of PFDGS to reflect ECM metabolism. Taken together, a high PFDGS score represented multiple adverse biological alterations, including a high level of immunoinflammatory activation, disordered peripheral blood myeloid cells and lymphocytes, and a high level of ECM hypermetabolism, which may contribute to the biological interpretation that patients with high PFDGS scores possessed a higher risk of IA onset and rupture.