Identification and validation of key genes mediating intracranial aneurysm rupture by weighted correlation network analysis.

BACKGROUND: Rupture of intracranial aneurysm (IA) is the leading cause of subarachnoid hemorrhage. However, there are few pharmacological therapies available for the prevention of IA rupture. Therefore, exploring the molecular mechanisms which underlie IA rupture and identifying the potential molecular targets for preventing the rupture of IA is of vital importance.
METHODS: We used the Gene Expression Omnibus (GEO) datasets GSE13353, GSE15629, and GSE54083 in our study. The 3 datasets were merged and normalized. Differentially expressed gene (DEG) screening and weighted correlation network analysis (WGCNA) were conducted. The co-expression patterns between ruptured IA samples and unruptured IA samples were compared. Then, the DEGs were mapped into the whole co-expression network of ruptured IA samples, and a DEG co-expression network was generated. Molecular Complex Detection (MCODE) (http://baderlab.org/Software/MCODE) was used to identify key genes based on the DEG co-expression network. Finally, key genes were validated using another GEO dataset (GSE122897), and their potential diagnostic values were shown using receiver operating characteristic (ROC) analysis.
RESULTS: In our study, 49 DEGs were screened while 8 and 6 gene modules were detected based on ruptured IA samples and unruptured IA samples, respectively. Pathways associated with inflammation and immune response were clustered in the salmon module of ruptured IA samples. The DEG co-expression network with 35 nodes and 168 edges was generated, and 14 key genes were identified based on this DEG co-expression network. The gene with the highest degree in the key gene cluster was CXCR4. All key genes were validated using GSE122897, and they all showed the potential diagnostic value in predicting IA rupture.
CONCLUSIONS: Using a weighted gene co-expression network approach, we identified 8 and 6 modules for ruptured IA and unruptured IA, respectively. After that, we identified the hub genes for each module and key genes based on the DEG co-expression network. All these key genes were validated by another GEO dataset and might serve as potential targets for pharmacological therapies and diagnostic markers in predicting IA rupture. Further studies are needed to elucidate the detailed molecular mechanisms and biological functions of these key genes which underlie the rupture of IA. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

Intracranial aneurysm (IA) is the irreversible dilation of cerebral arteries, and the rupture of IA is the primary cause of subarachnoid hemorrhage (1,2). The prevalence of unruptured IA is about 3–5% in the general population, and about 20–30% of these patients have multiple aneurysms (3-5). The main risk factors for IA rupture include hypertension, smoking, aneurysm size, and location. Currently, endovascular coiling and surgical clipping are the two main preventative methods for IA rupture (6). Although widely used, the overall morbidity and mortality associated with these two invasive procedures are not negligible, and pharmacological therapies might provide an alternative for those who have high risks for invasive procedures (7). However, there are few pharmacological therapies available for the prevention of IA rupture (6,7). Thus, exploring the molecular mechanisms underlying the event of IA rupture and identifying the potential molecular targets for preventing IA rupture are vitally important.

Previous studies have suggested that inflammation and immune response, extracellular matrix degradation, hemodynamic stress, etc. participate in the pathogenesis and rupture of IA (8-12). Although many studies have performed detailed experiments to elucidate the molecular mechanisms behind IA rupture, few studies have used a weighted correlation network analysis (WGCNA) approach based on high-throughput data (such as mRNA microarray and RNA-seq) to identify potential targets that might have a critical biological function in the mechanism of IA rupture.

The WGCNA algorithm was created by Zhang & Horvath in 2005 and was initially used to create a gene co-expression network from mRNA microarray data (13). WGCNA generates a scale-free gene co-expression network based on Pearson’s correlation matrix of genes and detects gene modules from the network generated. These gene modules can be used for further analysis, such as correlating modules with clinical traits, functional enrichment, and hub gene identification (13). Furthermore, WGCNA can be used to construct weighted correlation networks based on proteome, metabolome, and even microbiome data (14,15). The ‘WGCNA’ R package was available on the R official website (https://cran.r-project.org/) (16).

Using WGCNA based on two GEO datasets, Zheng et al. identified FOS, CCL2, COL4A, and CXCL5 as crucial genes mediating the pathogenesis of IA, and Liao et al. have conducted a similar analysis using another GEO dataset (17,18). Alongside crucial genes, Bo et al. identified several crucial miRNAs as mediators of IA pathogenesis (19). Additionally, Landry et al. combined WGCNA and population-specific gene expression analysis (PSEA) to analyze the microenvironment associated with ruptured IA (20). Although WGCNA has partially elucidated the pathogenesis of IA, the crucial genes mediating the rupture of IA require exploration.

In our study, we focused on key genes underlying the molecular mechanism of IA rupture. We included 3 mRNA microarray datasets from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Expression profiling data of unruptured IA and ruptured IA were used to construct the co-expression network. The gene modules were detected, and the genes in each module were subjected to functional enrichment analysis. Differentially expressed genes (DEGs) were mapped into the co-expression network and constructed into a DEG co-expression network. Key genes were identified based on this DEG co-expression network, and gene set enrichment analysis (GSEA) was conducted for them. Finally, the key genes were further validated by using another GEO dataset, and the potential diagnostic value of these key genes was determined. We present the following article in accordance with the MDAR reporting checklist (available at http://dx.doi.org/10.21037/atm-20-4083).

Methods

Medical ethics

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The raw datasets were available from the GEO database (http://www.ncbi.nlm.nih.gov/geo/; GSE13353, GSE15629, GSE54083 and GSE122897). In our study, neither human trials nor animal experiments were executed.

Data collection

The datasets GSE13353, GSE15629, and GSE54083 were selected from the GEO database. Series matrix files and data tables of these microarray platforms (GPL570, GPL6244, and GPL4133, respectively) were downloaded.

Data preprocessing and DEG screening

Microarray probes in the series matrix files were annotated with official gene symbols using the data table of corresponding platforms, and gene expression matrices were obtained. Then, these gene expression matrices were merged into one matrix, and the sva R package was used to normalize the expression data from 3 different datasets. The limma R package was used to conduct the DEG analysis based on the normalized gene expression matrix. The threshold for DEG analysis was set as |log2 (fold-change)| >0.5 and P<0.05.

Weighted co-expression network construction

The WGCNA R package was used to construct the co-expression networks of ruptured IA samples and unruptured IA samples. Genes with the top 25% variance were chosen, and 27 ruptured IA samples and 19 unruptured IA samples were involved in the two WGCNA analyses, respectively. The adjacency matrices which store the information of the whole co-expression network were created based on Pearson’s correlation matrices. Topological overlap measure (TOM) matrices were obtained from adjacency matrices to detect gene modules of the co-expression networks. Based on the TOM matrices, the average linkage hierarchical clustering method was used to construct clustering dendrograms with a minimum module size of 30. Finally, similar gene modules were merged, with a threshold of 0.25.

Functional annotation of gene modules

Genes in each module were subjected to gene ontology (GO) and Kyoto Encyclopedia of Genes Genomes (KEGG) pathway analysis to understand their biological function better. The Database for Annotation, Visualization, and Integrated Discovery (DAVID) v6.8 (http://david.abcc.ncifcrf.gov/) was used for analysis. The cut-off for functional annotation was set as count >2 and P<0.05. The co-expression patterns between ruptured IA samples and unruptured IA samples were compared based on the results of the functional annotations.

Hub gene identification and key gene mining

The intramodular connectivities of the genes were calculated using the WGCNA algorithm, and the genes with the highest intramodular connectivity were identified as hub genes for each module. Next, using Cytoscape v3.7.0 (https://cytoscape.org/), DEGs were mapped into the whole co-expression network of ruptured IA samples. The unruptured IA samples and DEG co-expression networks were obtained after removing isolated nodes and isolated node pairs. The Molecular Complex Detection (MCODE), a plugin in Cytoscape to detect densely connected regions in a given network, was used to generate the key gene cluster based on the DEG co-expression network of ruptured IA samples. Both the DEG co-expression network and hub gene cluster of ruptured IA samples were visualized.

GSEA and validation of key genes

The GSEA v4.0.3 software was used to conduct GSEA for each key gene. The c5.bp.v7.0.symbols.gmt downloaded from the GSEA official website (https://www.gsea-msigdb.org/gsea/downloads.jsp) was used as the reference gene sets. For each key gene, the IA samples were divided into two groups according to the median expression of the particular key gene to perform GSEA. Finally, another GEO dataset (GSE122897) containing RNA-seq data of ruptured IA samples and unruptured IA samples was chosen for validation of the key genes. Finally, key genes were used to conduct receiver operating characteristic (ROC) analysis to show the potential diagnostic value of these genes.

Statistical analysis

Data preprocessing, DEG screening and weighted co-expression network analysis were performed in R v3.6.2. Functional annotation analysis was performed by DAVID v6.8. Key genes were mined using MCODE in Cytoscape v3.7.0. GSEA was conducted by GSEA v4.0.3. The details of these bioinformatic analyses have been described in corresponding subsections. The potential diagnostic value of key genes was shown by ROC analysis using IBM SPSS 25.0. A P value <0.05 was considered statistically significant.

Results

Workflow

The workflow of our study is shown in . After merging 3 GEO datasets and batch normalization (data preprocessing), we conducted DEG screening and WGCNA analysis. Gene modules were detected based on the co-expression networks. GO, and KEGG pathway analyses were conducted to compare the co-expression patterns between ruptured IA and unruptured IA. After that we created a DEG co-expression network by mapping the screened DEGs into the whole co-expression network of a ruptured IA co-expression network. MCODE was used to mine key genes based on the DEG co-expression network, and ROC analysis was used to determine the potential clinical significance of key genes. Finally, single-gene GSEA and validation of key genes were performed.

Figure 1

Workflow of our study. DEG, differentially expressed gene; ROC, receiver operating characteristic; GSEA, gene set enrichment analysis.

DEG screening

We screened 49 DEGs, of which 28 genes were up-regulated, and 21 genes were down-regulated in ruptured IA samples. DEGs with top-10 fold-change are shown in . The expression patterns of DEGs between ruptured IA and unruptured IA are detailed in .

Table 1

DEGs with top-10 fold-change (ruptured IA/unruptured IA)

Gene symbol	Official full name	Log2(fold-change)
Up-regulated
PPBP	Pro-platelet basic protein	1.67
PF4	Platelet factor 4	1.19
S100A8	S100 calcium binding protein A8	1.09
FPR1	Formyl peptide receptor 1	1.01
C15orf48	Chromosome 15 open reading frame 48	0.91
NCF2	Neutrophil cytosolic factor 2	0.91
RGS18	Regulator of G protein signaling 18	0.90
UPP1	Uridine phosphorylase 1	0.89
CXCR4	C-X-C motif chemokine receptor 4	0.88
SLA	Src like adaptor	0.87
Down-regulated
NR1D2	Nuclear receptor subfamily 1 group D member 2	−1.05
ATP1A2	ATPase Na+/K+ transporting subunit alpha 2	−0.95
FMO2	Flavin containing dimethylaniline monoxygenase 2	−0.89
SLC6A1	Solute carrier family 6 member 1	−0.77
CYP4X1	Cytochrome P450 family 4 subfamily X member 1	−0.76
NTRK3	Neurotrophic receptor tyrosine kinase 3	−0.73
PER3	Period circadian regulator 3	−0.66
ZIC1	Zic family member 1	−0.63
OPCML	Opioid binding protein/cell adhesion molecule like	−0.63
ZMAT1	Zinc finger matrin-type 1	−0.61

DEG, differentially expressed gene; IA, intracranial aneurysm.

Figure 2

DEG expression patterns. (A) Heatmap for DEGs. Up-regulated genes were in red color while down-regulated genes were in green color. (B) Volcano plot for DEGs. Red dots represent up-regulated genes while green dots represent down-regulated genes. DEG, differentially, expressed gene; FC, fold change.

Weighted gene co-expression network construction

For the ruptured IA samples, 3 outliers (GSM337109, GSM337119, and GSM391293) were detected (), and for the unruptured IA samples, 1 outlier (GSM337112) was detected (Figure S1A). After removing the outliers, 24 ruptured IA samples and 18 unruptured IA samples were included in further analysis. The soft-threshold power β=16 and 18 were chosen for the ruptured IA samples and unruptured IA samples, respectively ( & Figure S1B,C,D,E). Finally, we detected 8 gene modules for ruptured IA samples and 6 gene modules for unruptured IA samples (, Figure S1F).

Figure 3

Construction of weighted gene co-expression network of ruptured IA samples. (A) Sample clustering detected 3 outliers (GSM337109, GSM337119 and GSM 391293). (B,C) The cut-off was set to be 0.8 and β=16 was chosen to be the soft-threshold power. (D,E) The histogram and the linear fitting plot showed that the co-expression network we constructed met the requirement of scale-free topology. (F) 8 modules were detected after merging modules generated by dynamic tree cut. IA, intracranial aneurysm.

Co-expression pattern comparison

The results of GO-BP and KEGG analysis are shown in Table S1. The salmon module of ruptured IA samples was mainly associated with inflammation and immune response. The GO-BP terms and KEGG pathways with top-10 count numbers for ruptured IA samples are shown in and . However, these inflammation and immune-related pathways were scattered in different modules of unruptured IA samples.

Figure 4

GO-BP and KEGG pathway for salmon module of ruptured IA. GO, gene ontology; BP, biological process; KEGG, Kyoto encyclopedia of genes and genomes; IA, intracranial aneurysm.

Table 2

GO-BP terms and KEGG pathways with top-10 count number for ruptured IA

ID	Terms	Count	-logP
GO-BP
GO:0007165	Signal transduction	13	1.94
GO:0006954	Inflammatory response	11	4.77
GO:0006955	Immune response	10	3.64
GO:0045087	Innate immune response	9	2.89
GO:0042742	Defense response to bacterium	7	4.09
GO:0007155	Cell adhesion	7	1.58
GO:0002576	Platelet degranulation	6	3.81
GO:0051092	Positive regulation of NF-kappaB transcription factor activity	6	3.29
GO:0032496	Response to lipopolysaccharide	6	2.88
GO:0030593	Neutrophil chemotaxis	5	3.50
KEGG pathway
hsa04060	Cytokine-cytokine receptor interaction	8	2.93
hsa04062	Chemokine signaling pathway	5	1.43

GO, gene ontology; BP, biological process; KEGG, Kyoto encyclopedia of genes and genomes; IA, intracranial aneurysm; NF, nuclear factor.

Hub gene and key gene mining

The WGCNA algorithm calculated the intramodular connectivity of genes in each module, and the gene with the highest intramodular connectivity was considered to be the hub gene of that module (Table S2). Thus, vanin 2 (VNN2), a pro-inflammatory gene, was deemed the hub gene for the salmon module of the ruptured IA group.

A DEG co-expression network with 35 nodes and 168 edges was generated by mapping DEGs into the whole co-expression network of DEGs (). MCODE obtained key genes based on this DEG co-expression network ( & ), of which 12 were up-regulated, and 2 were down-regulated. The gene with the highest degree in the key gene cluster was C-X-C motif chemokine receptor 4 (CXCR4).

Figure 5

DEG-co-expression network and key gene cluster. (A) Up-regulated genes were in red while down-regulated genes were in green. (B) Key gene cluster detected by MCODE. DEG, differentially expressed genes; MCODE, molecular complex detection.

Table 3

Key genes mediating rupture of IA

Gene ID	Gene symbol	Official full name
84419	C15orf48	Chromosome 15 open reading frame 48
366	AQP9	Aquaporin 9
6503	SLA	Src like adaptor
4354	MPP1	Membrane palmitoylated protein 1
23024	PDZRN3	PDZ domain containing ring finger 3
284	ANGPT1	Angiopoietin 1
2357	FPR1	Formyl peptide receptor 1
290	ANPEP	Alanyl aminopeptidase, membrane
5552	SRGN	Serglycin
4688	NCF2	Neutrophil cytosolic factor 2
128346	C1orf162	Chromosome 1 open reading frame 162
10312	TCIRG1	T cell immune regulator 1, ATPase H+ transporting V0 subunit a3
7852	CXCR4	C-X-C motif chemokine receptor 4
2268	FGR	FGR proto-oncogene, Src family tyrosine kinase

IA, intracranial aneurysm.

In the GSEA results for each key gene, we selected and visualized the significantly enriched gene sets with the highest normalized enrichment score (NES) value in each analysis (). Most of the gene sets that we selected were associated with the inflammatory or immune responses, except for 3 gene sets, of which their related genes were MPP1, FPR1, and TCIRG1.

Figure 6

Single gene GSEA for key genes. Most GSEA results of key genes were associated with inflammation/response except 3 key genes (MPP1, FPR1 and TCIRG1). GSEA, gene set enrichment analysis.

These key genes were validated using GSE122891, and the expression level of all key genes showed significant differences between the ruptured IA samples and unruptured IA samples ( & Table S3).

Figure 7

Validation of key genes using dataset GSE122897. The expression levels of key genes were compared using unpaired t-test and the expression level of all these genes showed significant difference between ruptured IA samples and unruptured IA samples. IA, intracranial aneurysm. *, P<0.05; **, P<0.01; ***, P<0.001.

Finally, the potential clinical values of key genes were shown by ROC analysis. All the genes showed potential diagnostic values (P<0.05) (Table S4). ROC curves of key genes with the top-6 degree were visualized ().

Figure 8

ROC curves for key genes with top-6 degree in key gene cluster. ROC, receiver operating characteristic; AUC, area under curve.

Discussion

In our research, 49 DEGs were screened, among these, 28 genes were up-regulated, and 21 genes were down-regulated in the ruptured IA samples ( & ). The weighted gene co-expression networks for both ruptured IA and unruptured IA were constructed. Based on the co-expression networks, 8 and 6 gene modules were detected for ruptured IA samples and unruptured IA samples, respectively (, Figure S1). Genes in each module were subjected to GO and KEGG pathway enrichment analysis and co-expression patterns between ruptured IA samples. Pathways related to inflammation and immune response were clustered in the salmon module of ruptured IA samples, yet scattered in different modules of unruptured IA samples ( & ), and the hub gene of the salmon module of ruptured IA samples was VNN2, which is a pro-inflammatory gene (21). Then, a DEG co-expression network with 35 nodes and 168 edges was generated by mapping DEGs into the whole co-expression network of ruptured IA samples (). MCODE obtained key genes based on the DEG co-expression network ( & ). Single gene GSEA showed that the gene sets most enriched in groups divided by the expression level of key genes were associated with inflammation and immune response (). Finally, all these key genes were validated by another GEO dataset ( & Table S3), and they also showed potential diagnostic value ( & Table S4).

Kleinloog summarized the DEGs that were screened by their research and other previous studies. Almost all of the key genes in our study (except PDZRN3 and ANGPT1) were following the DEGs reported by the previous studies (22-26). However, some DEGs were not overlapping with DEGs in previous studies. The difference might be the result of the datasets selected, different thresholds, and different algorithms for the DEG screening.

Previous and recent studies have indicated that inflammatory and immune response played a critical role in the rupture of IA. Using high throughput data (mRNA microarray and RNA-seq), Kleinloog, Pera, and Nakaoka et al. discovered that inflammation and immune response participated in the rupture of IA (22,24,26). Inflammatory and immune cells such as macrophages, T, and B lymphocytes were also found to be potential mediators of IA rupture. Kataoka et al. observed significantly increased macrophage infiltration in ruptured human IA (27). In lymphocyte deficient and wild-type mice models, Sawyer et al. discovered that both the formation of IA and rupture of IA were significantly fewer in lymphocyte deficient mice compared with wild-type mice (28). Also, pro-inflammatory factors have critical roles in the progression of IA rupture. A well-known pro-inflammatory cytokine, TNF-α, can cause the weakening of arterial walls through downstream pathways associated with endothelial dysfunction, vascular smooth muscle phenotypic modulation, activating pro-inflammatory cells, and other mechanisms (29,30). Another inflammatory mediator, NF-κB, participates in both the formation and rupture of IA (31,32). Our finding that the pathways associated with inflammation and immune response were clustered in the salmon module of ruptured IA samples, and our GSEA results that inflammation and immune-related gene sets were up-regulated, were consistent with these previous and recent studies.

VNN2, also named GPI-80, was the hub gene in the salmon module of ruptured IA samples. It was first discovered on human neutrophils, and it regulates neutrophil adhesion and transendothelial migration (33,34). VNN2 was also found on a subpopulation of monocytes with a superior ability of reactive oxygen species production and phagocytosis (35). Being that neutrophils and monocytes are involved in inflammation and participate in the formation and rupture of IA (12,36,37), VNN2 might have potential roles in mediating IA rupture.

One of the widely studied chemokine receptors, CXCR4, was the gene with the highest degree in the key gene cluster. CXCR4 was firstly described as a G protein coupled-receptor that mediates the fusion of HIV-1 with host cells, and CXCL12 and SDF-1 were the specific ligands for CXCR4 (38-40). The signaling pathway of CXCR4 plays various roles in essential biological processes such as vascular development, cardiogenesis, hematopoiesis, and inflammation, among others (41).

Previous studies have described the association between CXCR4 and pathophysiology of the aneurysm. In mice models with carotid or IA, Hoh et al. found that expression of CXCR4 is increased in circulating progenitor cells, which may later differentiate into inflammatory cells (42). Besides, in the abdominal aortic aneurysm (AAA), research has shown that CXCR4 mediated aneurysmal progression, possibly through the inflammatory response. Michineau et al. found that SDF-1α and CXCR4 were up-regulated in CaCl2 induced AAA mice and that the block of CXCR4 by AMD3100 inhibited AAA formation and progression. They also observed decreased infiltration of macrophages along with decreased levels of pro-inflammatory factors such as MCP-1, IL-1β, and TNF-α (43). Another group discovered that CXCR4 was highly expressed in T/B lymphocytes and macrophages in AAA, also suggesting that CXCR4 promotes AAA formation and progression through the inflammatory response (44). Kodali et al. also found that metalloproteinase-2, implicated in aneurysm formation and progression, could be induced by SDF-1 via binding to CXCR4 (45). Our results that CXCR4 was the up-regulated key gene with the highest degree and inflammatory and immune response-related gene sets were up-regulated in the CXCR4 high group were consistent with these previous studies. Furthermore, CXCR4, as well as other key genes, also showed potential diagnostic values and might serve as biomarkers to predict IA rupture.

In our study, we constructed co-expression networks for ruptured IA and unruptured IA samples using 3 GEO datasets by WGCNA, detected gene modules, and identified key genes mediating the rupture of IA for the first time. Also, we conducted a GSEA analysis for the key genes and validated the key genes we identified using another GEO dataset. However, our study had some limitations. Due to the lack of detailed information about GEO datasets, it is difficult for us to associate the gene modules and key genes with clinical traits.

Our study revealed that the co-expression pattern between ruptured IA and unruptured IA was different. These findings implicated that the inflammatory and immune response might play important roles in IA rupture. The key genes that we identified in this study might have crucial biological functions in mediating rupture of IA.

Conclusions

Using a weighted gene co-expression network approach, we identified 8 and 6 modules for ruptured IA and unruptured IA, respectively. Next, we identified the hub genes for each module and the key genes based on the DEG co-expression network. These key genes were validated by another GEO dataset (GSE122897) and might serve as potential targets for pharmacological therapies and diagnostic markers for predicting IA rupture. Further studies are needed to elucidate the detailed molecular mechanisms and biological functions of these key genes that underlie the rupture of IA.

The article’s supplementary files as