The workflow chart showing the study design is presented in Fig. 1. Five datasets of IA were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), including four microarray datasets (GSE75436, GSE54083, GSE13353, GSE15629) and one RNA sequencing (RNA-seq) dataset (GSE122897). Details of the selected datasets are shown in Supplementary Table 1. Raw data from the four microarrays were merged. Batch effects were removed by package “SVA” [20], and data normalization was implemented by package “limma” in R [21] (Supplementary Fig. 1). Corresponding clinical data were obtained from the appropriate metadata within the GEO database. The merged microarray data were used as the training set, and the RNA-seq data used as the validation set.

Fig. 1.

Workflow chart. IA, intracranial aneurysm; DEG, differentially expressed gene; WGCNA, weighted gene co-expression network analysis.

Differentially expressed genes (DEGs) were identified using the empirical Bayesian approach of the “limma” package [21]. Genes with p 0.05 and $|$log2 fold change (FC)$|$ $>$0.5 between IA samples and control artery samples were considered to be DEG. The “pheatmap” package was used to draw heatmaps, and the “ggplot2” package used to draw volcano plots. The overlap of DEGs between the microarray and the RNA-seq datasets was shown using the “VennDiagram” package.

Gene set enrichment analysis (GSEA) is a method to assess the distribution trend of genes in a predefined set [22]. In the present study, GSEA was used to assess the contribution of different cell death pathways to IA formation. The apoptosis, necroptosis, autophagy and ferroptosis pathways were downloaded from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp/), and the pyroptosis gene sets were downloaded from the Molecular Signatures Database (MSigDB) (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). GSEA was performed using the “ClusterProfiler” package and visualized using the “enrichplot” package. The threshold was set to an adjusted p (p.adj) 0.05. Gene set variation analysis (GSVA) was performed using the “GSVA” package to quantify the expression levels of various cell death pathways in all samples [23]. A p 0.05 was considered to indicate statistical significance between the IA and control groups. We also used the “Pathview” R package to visualize the relative expression levels of key genes in the KEGG pathways of apoptosis, necroptosis, autophagy and ferroptosis in IA samples [24]. In addition, the major positive and negative regulators of apoptosis, necroptosis, autophagy, ferroptosis and pyroptosis were obtained from the previous literature (Supplementary Table 2) [13, 14, 19]. The different expression levels of these regulators between the IA and control samples were visualized by heatmaps.

Enrichment analysis suggested that multiple cell death pathways were enriched in the IA samples. Therefore, we constructed an integrative classification of regulated cell death. The sets of apoptosis, necroptosis and pyroptosis genes were downloaded from MSigDB, the autophagy-related genes were downloaded from the Human Autophagy Database (HADb) (http://www.autophagy.lu/), and the ferroptosis-related genes were downloaded from FerrDb V2 (http://www.zhounan.org/ferrdb/current/). The various types of cell death-related genes were intersected with DEGs separately to obtain the cell death signature genes.

To further identify different IA patterns related to cell death signature genes, an unsupervised consensus clustering algorithm was applied with the k-means method using the “ConsensusClusterPlus” package [25]. Clustering was iterated 1000 times to ensure the stability of the classification. The optimal number of clusters was determined by the relative changes in the area under the cumulative distribution function (CDF) curves of the consensus score and the consensus heatmap. The novel cell death classification was visualized and verified using a principal component analysis (PCA) with the “factoextra” package. Single sample GSEA (ssGSEA) was applied to quantify the enrichment level of the set of cell death signature genes for each sample using the “GSVA” package. ssGSEA scores were defined as cell death scores.

DEGs were screened using the “limma” package in two cell death clusters. p 0.05 and $|$log2 FC$|$ $>$1 were defined as the threshold for differential expression of genes between two clusters. The results of differential gene expression were displayed using volcano graphs and heatmaps. The overlapping DEGs identified in the microarray and RNA-seq datasets were shown in Venn diagrams.

The “ClusterProfiler” R package was used to perform Gene Ontology (GO) functional and KEGG pathway enrichment analysis of the newly identified DEGs [26]. Significantly enriched functions and pathways were selected with an adjusted p-value (p.adj) 0.05 and false discovery rate (FDR) 0.2. To comprehensively explore differentially expressed levels of various signaling pathways in the two cell death patterns, the pathways for cell growth and death, signal transduction, signaling molecules and interaction, and immune system were downloaded from the KEGG database. The GSVA was implemented and the GSVA scores of different signaling pathways between the two cell death clusters were compared using the “limma” tool. A p.adj 0.05 was considered to indicate statistical significance. Furthermore, the differentially expressed levels of the main regulators of various cell death types between two clusters were compared by the “limma” package.

The algorithms “MCPcounter” [27] and “xCell” [28] were used to examine the immunocyte infiltration of IAs. The immune cell composition of patients with different cell death clusters was compared and shown in boxplots. Correlations between cell death scores and the infiltration levels of different immune cell types were analyzed using Spearman correlation analysis.

Weighted gene co-expression network analysis (WGCNA) was performed using the “WGCNA” R package to determine the gene modules linked to the different cell death patterns [29]. The WGCNA network was first constructed using a scale-free topology network with an optimal soft threshold. The relationship between the gene modules and two cell death clusters was then evaluated by Spearman correlation analysis and shown in heatmaps. The genes in the module with the strongest correlation with cluster A were imported into the STRING database (https://cn.string-db.org/). The protein-protein interaction (PPI) network was obtained with high confidence (minimum required interaction score = 0.700) and visualized using Cytoscape software (version 3.9.1, https://cytoscape.org/, GitHub Inc, San Francisco, CA, USA). Finally, the CytoHubba plugin was used to determine the hub proteins in the PPI network via the degree algorithm.

The merged microarray dataset was used as the training set, and the RNA-seq dataset was used as the validation set. The DEGs identified between two cell death clusters were used for screening. First, the least absolute shrinkage and selection operator (LASSO) technique and random forest (RF) algorithm were used to eliminate redundant variables. The optimal parameter with the lowest error rate and the best stability tree number was determined by calculating the error rates for 1 to 1000 trees. The Gini coefficient method was applied to evaluate the importance of the variables through RF. Univariate logistic regression was then performed to further select valuable predictive genes for IA rupture. The generalized linear model was finally used to develop a prediction model of IA rupture. The predictive performance of the generalized linear model was measured by the concordance index (C index) and receiver operating characteristic (ROC) curves [30]. Calibration was evaluated using a calibration plot [31]. The goodness of fit of the model was assessed by the Hosmer-Lemeshow (HL) test. LASSO regression, logistic regression, and generalized linear model were conducted using the “glmnet” R package. The random forest was carried out using the “random Forest” package in R. ROC curves were performed using the “pROC” package.

All of the statistical analysis was performed in R (version 4.2.1, https://www.r-project.org, the Comprehensive R Archive Network, Vienna, Austria). Continuous data were expressed as means $\pm$ standard deviations (SDs), or quartiles (medians [25th percentiles, 75th percentiles]), while categorical data were summarized as frequencies and percentages. Student’s t test or the Wilcoxon test were used for pairwise comparison between two groups. The statistical significance in all of the analyses was assumed to be p 0.05.

The merged microarray dataset included 61 IA and 30 control samples. The RNA-seq dataset included 44 IA and 16 control samples. The microarray cohort contained 27 ruptured IAs and 19 unruptured IAs, with the rupture status of the remaining 15 IA samples being unknown. The RNA-seq cohort contained 21 ruptured IAs, 21 unruptured IAs, and two IAs with unclear rupture status. A total of 1849 DEGs were identified in the microarray dataset, including 975 upregulated genes and 874 downregulated genes. In the RNA-seq dataset, 1486 upregulated DEGs and 1506 downregulated DEGs were identified. Volcano plots and heatmaps were used to visualize the expression of DEGs in all cases (Fig. 2A,B). The Venn diagrams show the intersection of DEGs between the microarray and RNA-seq datasets. These yielded 350 upregulated DEGs and 257 downregulated DEGs (Fig. 2C).

Fig. 2.

Identification of differentially expressed genes (DEGs) and analysis of cell death enrichment. (A) Volcano plot and heatmap of DEGs identified in the microarray set. (B) Volcano plot and heatmap of DEGs identified in the RNA-seq set. (C) Venn diagram of DEGs. (D,E) Gene set variation analysis (GSVA) of cell death pathways. (F,G) Gene set enrichment analysis (GSEA) of cell death pathways.

GSEA and GSVA were performed to explore the role of various types of RCD in IA formation. GSVA results in the microarray cohort showed that apoptosis, necroptosis, and ferroptosis expression levels were significantly higher in the IA samples than in controls (Fig. 2D). In the RNA-seq dataset, apoptosis, necroptosis, and pyroptosis showed higher expression in IA than in control samples (Fig. 2E). In the microarray dataset, GSEA revealed that IAs were significantly enriched in apoptosis, ferroptosis and pyroptosis-related pathways (Fig. 2F). In the RNA-seq set, DEGs were markedly enriched in apoptosis, necroptosis and pyroptosis (Fig. 2G). The relative expression levels of key genes in the apoptosis, necroptosis, autophagy, and ferroptosis pathways were visualized using pathview plots (Supplementary Fig. 2). The expression levels of multiple positive regulators of cell death were higher in IAs than in controls, while several negative regulators were downregulated in IAs (Fig. 3A,B).

Fig. 3.

Expression of cell death-related regulators. (A,B) Heatmaps visualizing the differential expression levels of cell death-related regulators between IA and control tissues. (C,D) Heatmaps visualizing the differential expression levels of cell death-regulators between the two cell death clusters. *p 0.05, **p 0.01, ***p 0.001.

Functional enrichment analysis showed that multiple types of cell death were involved in IA. Therefore, we further investigated the cell death signature genes associated with IA. Firstly, we downloaded 161 apoptosis-related genes, 52 necroptosis-related genes, and 23 pyroptosis-related genes from the MSigDB database, as well as 222 autophagy-related genes from the HADb database and 198 ferroptosis-related genes from the FerrDb V2 database. Next, we obtained the overlapping genes by intersecting various cell death-related genes and the DEGs between IAs and controls (Fig. 4A). In total, we identified 15 signature genes related to apoptosis, 5 to necroptosis, 11 to autophagy, 8 to ferroptosis, and 2 to pyroptosis (Fig. 4A). After the removal of duplicates, a total of 33 cell death signature genes were identified. All showed increased expression in IA samples compared to controls, except for ITPR1, EGFR, LIFR and CPEB1 (Fig. 4B,C). Spearman correlation analysis also showed the expression of most of these cell death signature genes was positively correlated with each other (Supplementary Fig. 3). Signature genes were uploaded to the STRING database and the PPI network was visualized using Cytoscape software (Fig. 4D). The top 10 hub genes identified by CytoHubba plugin using the degree algorithm were EGFR, CXCR4, CTSB, FASLG, TIMP1, GZMB, NLRP3, CCR2, CD69 and CD2 (Fig. 4E).

Fig. 4.

Identification of cell death signature genes. (A) Venn diagram of various cell death-related genes and DEGs. Table showing cell death signature genes. (B,C) Boxplots of cell death signature genes in IAs compared to control tissues. (D) The PPI network of cell death signature genes. Orange represents upregulated DEGs, while blue represents downregulated DEGs. (E) Construction of PPI network for the top 10 cell death signature genes; the redder the color, the higher the rank. ns p $\mathrm{\ge }$ 0.05, *p 0.05, **p 0.01, ***p 0.001.

Consensus clustering analysis was used to determine the potential integrative cell death classification of IAs. According to the relative change in the area under the CDF curve and the consensus heatmap, the optimal number of clusters was determined to be two. All patients in the microarray dataset were divided into two cell death clusters, with no appreciable increase detected in the area under the CDF curve (Fig. 5A–C). A similar result was also obtained with the RNA-seq dataset (Fig. 5G–I). PCA showed the cell death signature could separate IAs into two distinguishable patterns (Fig. 5D,J). Based on the cell death signature ssGSEA scores, cluster A could be defined as a cell death-high subtype, and cluster B as a cell death-low subtype (Fig. 5E,K). Furthermore, the cell death-based classification was associated with different clinical phenotypes: cluster A had a high risk of rupture, but cluster B had a relatively lower rate of rupture (Fig. 5F,L). A consensus clustering analysis was also performed to identify potential IA subtypes based on different types of cell death. Apoptosis-, necroptosis-, autophagy-, ferroptosis- and pyroptosis-based clustering separated the IAs into two similar clusters. Overall, the Sankey diagrams indicated the cell death signature-based classification provided a better reflection of the clinical phenotypes than the clusters established for each independent type of cell death (Supplementary Fig. 4).

Fig. 5.

Novel cell death classification of intracranial aneurysms (IAs) determined from the microarray set (A–F) and from the RNA-seq set (G–I). (A,G) Cumulative distribution function (CDF) curves of the consensus score (k = 2–9). (B,H) Relative change in the area under the CDF curve (k = 2–9). (C,I) Consensus clustering matrix for k = 2, which was the optimal number of clusters. (D,J) Principal component analysis (PCA) according to cell death clusters. (E,K) Boxplots of single sample GSEA (ssGSEA) scores of IA samples in cluster A compared to cluster B. (F,L) Distribution of different rupture status groups within cell death clusters.

Given that patients in cluster A had a high risk of rupture, we further investigated DEGs and functional pathways between the two cell death clusters so as to better understand the potential mechanism of IA progression. DEGs between the two clusters are shown in Fig. 6A,B and Fig. 6E,F. GO and KEGG enrichment analysis revealed the DEGs were significantly enriched in immune-related terms (e.g., lymphocyte activation, mononuclear cell proliferation, macrophage activation, and mast cell activation). Furthermore, inflammatory pathways (e.g., cytokine activity and chemokine signaling pathway), ECM degradation-related processes (e.g., extracellular matrix disassembly and collagen catabolic process) were activated, but smooth muscle contraction and muscle tissue development were significantly inhibited (z-score 0) in cluster A (Supplementary Fig. 5C–H).

Fig. 6.

Subgroup analysis between the two cell death clusters in the microarray (A–D) and RNA-seq sets (E–H). Volcano plots (A,E) and heatmaps (B,F) of DEGs between the two cell death clusters. (C,G) Heatmaps visualizing the GSVA results of cell death, signal transduction, and immune system-related pathways between the different cell death clusters. (D,H) Spearman correlation analysis of the GSVA scores for these pathways in IAs. *p 0.05, **p 0.01, ***p 0.001.

GSVA revealed that patients in cluster A showed upregulation of cell growth and death-related pathways including apoptosis, ferroptosis, necroptosis and pyroptosis, as well as signal transduction pathways such as MAPK, VEGF, JAK-STAT, NF-$\kappa$B, TNF and HIF-1. Cluster B showed marked upregulation of Hedgehog, TGF-$\beta$ and Hippo signaling pathways. Samples in cluster A were also significantly enriched in various immune system-related pathways (Fig. 6C,G). Spearman correlation analysis showed that activation of apoptosis, necroptosis, ferroptosis and pyroptosis were strongly and positively correlated with the activation of VEGF, JAK-STAT, NF-$\kappa$B, TNF, Toll-like receptor, NOD-like receptor and chemokine signaling pathways. Furthermore, the expression of mTOR signaling pathways was significantly and positively correlated with the activation of autophagy (Fig. 6D,H).

Cluster A also showed overexpression of multiple apoptosis-positive regulators (e.g., CASPs, BAK1, BAX, BCL2L11, PMAIP1, BBC3, BID), necroptosis-positive regulators (RIPK3 and MLKL), autophagy-positive regulator (ULK2), ferroptosis-positive regulators (e.g., TFRC, ACSL4, LPCAT3, DPP4, NCOA4 and CARS1) and pyroptosis-positive regulators (e.g., CASP1, CASP4, CASP5 and GSDMD), but low expression of apoptosis-negative regulators (e.g., BCL2, BCL22L1, MCL1 and BCL2L2), necroptosis-negative regulator (PPM1B) and ferroptosis-negative regulators (e.g., HSPB1 and NFS1) (Fig. 3C,D). The simultaneous upregulation of multiple positive regulators and downregulation of negative regulators indicated that various cell death processes were activated in cluster A.

Differences in immune cell infiltration between the two cell death patterns was also investigated. In general, IA tissues in cluster A had higher immune scores and microenvironment scores than those in cluster B. MCPcounter analysis revealed that IA samples in cluster A had significantly higher levels of T cells, cytotoxic lymphocytes, B lineage cells, NK cells, monocytic lineage cells, and myeloid dendritic cells. xCell analysis also showed that IAs in cluster A had more infiltration of dendritic cells (DCs), mast cells, monocytes, macrophages and B cells (Fig. 7A). In addition, the infiltration levels of macrophages, M1 macrophages, mast cells, T cells and B cells were positively correlated with the cell death scores of IA samples (Fig. 7B). Similar results were observed with the RNA-seq dataset (Fig. 7C,D).

Fig. 7.

Immune characteristics of the two cell death clusters. (A,C) Immune cell content histogram. (B,D) Spearman correlation analysis between cell death scores and the infiltration levels of immune cells. ns p $\mathrm{\ge }$ 0.05, *p 0.05, **p 0.01, ***p 0.001.

We further analyzed correlations between the GSVA scores of different RCD and the degree of immune cell infiltration in all IA samples (Supplementary Fig. 6). This analysis showed that activation of apoptosis, necroptosis, ferroptosis and pyroptosis in IAs was positively correlated with the degree of infiltration by monocytes, macrophages, T cells and B cells. Interestingly, activation of the above types of RCD was significantly and positively correlated with the expression of M1 macrophages. However, the activation of M2 macrophages was positively correlated with activation of autophagy, indicating that autophagy plays an opposite role to other types of RCD in the regulation of macrophage phenotype switching.

The optimal soft threshold values in the microarray and RNA-seq cohorts were 3 and 6, respectively (Fig. 8A and Supplementary Fig. 7A). Fig. 8B and Supplementary Fig. 7B show the results of gene clustering and recognition of co-expression modules. In the microarray set, the turquoise module correlated significantly with cluster A (r = 0.86, p 0.001), while the blue module was significantly related to cluster B (r = 0.45, p 0.001) (Fig. 8C). In the RNA-seq cohort, the green/yellow module was significantly associated with cluster A (r = 0.74, p 0.001), and the cyan module to cluster B (r = 0.85, p 0.001) (Supplementary Fig. 7C). The genes within the modules were significantly correlated with the traits (cluster A and cluster B) (Fig. 8D,E and Supplementary Fig. 7D,E). The genes in the turquoise module were uploaded to the STRING database and the degree scores of each protein calculated by Cytoscape using the CytoHubba plugin. The PPI network found the core proteins related to cluster A were CD4, TYROBP, ITGB2, CD8A, SYK, IL10, TLR2, CD86, CCR5, CDF1R, LYN, LCP2, CD44, VAV1 and CD28 (Fig. 8F). A similar PPI network was obtained from the green/yellow module in the RNA-seq cohort (Supplementary Fig. 7F).

Fig. 8.

Weighted gene co-expression network analysis (WGCNA) in the microarray cohort. (A) Scale-free fit index and mean connectivity for different soft powers (optimal soft threshold $\beta$ = 3). (B) Gene clustering and recognition of co-expression modules. (C) Heatmap of the correlations between gene modules and cell death clusters. (D,E) Correlation between traits (cluster A and cluster B) and genes in the modules. (F) Construction of the PPI network of hub genes in the turquoise module: the redder the color, the higher the rank.

Since IAs in cell death cluster A had a relatively higher risk of rupture, DEGs between cell death cluster A and B were used to identify predictive variables for rupture. Following the removal of samples with unclear rupture status, 46 cases were obtained from the microarray dataset (as the training cohort) and 42 from the RNA-seq dataset (as the validation cohort). LASSO regression in the training cohort retained 23 genes (Fig. 9A,B). In addition, the top 20 genes according to the mean decrease Gini by RF were also selected as candidate variables (Fig. 9C,D). Seven genes were thus identified by intersecting the candidate genes obtained by the LASSO and RF algorithms (Fig. 9E). Subsequently, univariate regression identified ANPEP, STK10, DSC2, HPSE, IGFBP1, and S100A8 as risk factors for IA rupture, each with an odds ratio (OR) $>$1 and a p 0.05 (Fig. 9F). ROC curve analysis revealed the efficacy of these 6 genes for independently predicting IA rupture (Fig. 9G). Boxplots showed that all of them were upregulated in ruptured IAs (Fig. 9H). Finally, a prediction model was constructed based on the 6 characteristic genes and using a multifactor, generalized linear model: Y = (–37.5607) + (1.3321) $\times$ ANPEP + (–0.694) $\times$ STK10 + (0.6081) $\times$ DSC2 + (1.93) $\times$ HPSE + (1.4639) $\times$ IGFBP1 + (0.7208) $\times$ S100A8. The model demonstrated good accuracy in estimating the risk of rupture, with a C index of 0.92 (95% CI, 0.85–1.00) and an area under curve (AUC) of 0.924 in the training cohort (Fig. 9I). The p value of the HL test was 0.841 ($>$0.05) and the calibration curve was close to 45${}^{\circ }$, indicating a well-calibrated model (Fig. 9J). In the validation cohort, the predictive model displayed a C index of 0.85 (95% CI, 0.74–0.97) and an AUC of 0.855 for the prediction of IA rupture. The calibration plot also showed accurate predictive ability in the validation cohort (Supplementary Fig. 8).

Fig. 9.

Construction of the prediction model for IA rupture in the training cohort. (A) Screening of candidate genes using the LASSO algorithm. (B) Cross-validation for tuning parameter screening in the LASSO regression model. (C) Correlation plot between the number of trees in random forest (RF) and the model error. (D) The mean decrease in Gini index of the top 20 genes calculated in the RF classifier. (E) Venn diagram of candidate genes obtained from LASSO and RF. (F) Univariate logistic regression analysis of the seven candidate genes, which yielded six genes linked to rupture. (G) ROC curve and AUC values of the six risk genes for the prediction of rupture. (H) Boxplots of six risk genes in ruptured IAs compared to unruptured IAs. (I) ROC curve and AUC values of the prediction model in rupture prediction. (J) Calibration plot of the prediction model. *p 0.05, **p 0.01, ***p 0.001.