Integrated analysis of microarray data to identify the genes critical for the rupture of intracranial aneurysm.

Intracranial aneurysm (IA) is a localized dilation of the blood vessel. The present study was designed to explore the mechanisms of rupture of IA. GSE13353 (including 11 ruptured and 8 unruptured IA samples) and GSE15629 (including 8 ruptured and 6 unruptured IA samples) were downloaded from the Gene Expression Omnibus database. The differentially expressed genes (DEGs) identified using limma and MetaDE packages were merged, and a protein-protein interaction (PPI) network analysis was performed using Cytoscape software. Pathway enrichment analysis was performed for the nodes of the PPI network using the fisher algorithm. The 100 most prominent genes in the network were designated candidate genes and a hierarchical clustering analysis was performed. The tune.svm function of e1071 package was used to construct a support vector machine (SVM) classifier, and the Candidate Cancer Gene Database was applied to analyze the characterization of gene-associated cancer. Furthermore, the genes involved in the SVM classifier were assessed via principal component analysis (PCA). In the ruptured samples, 1,292 DEGs and 1,029 DEGs separately were identified by limma and MetaDE packages. The 100 most prominent genes in the network included fibronectin 1 (FN1), amyloid β (A4) precursor protein (APP), nuclear RNA export factor 1 (NXF1) and signal transducer and activator of transcription 3 (STAT3). Pathway enrichment analysis identified that toll-like receptor 3 (TLR3) was enriched in the Toll-like receptor signaling pathway. A total of 15 genes (including FN1) were used to construct the SVM classifier. NXF1 was identified to be associated with Nervous System Cancer. PCA revealed that APP, NXF1 and STAT3 were the 3 principal components. TLR3, FN1, APP, NXF1 and STAT3 may affect the rupture of IA.

Introduction

As a cerebrovascular disorder, intracranial aneurysm (IA; also designated brain or cerebral aneurysm) is a ballooning or localized dilation of the blood vessel induced by weakness of the artery wall (1). IA cases are divided according to size into small (diameter, <15 mm) and large aneurysms [including large (15–25 mm), giant (25–50 mm) and super-giant (>50 mm) aneurysms] (2). Based on the shape, IA may be classified into saccular aneurysms, fusiform aneurysms and microaneurysms (2). IA may not only be a result of genetic conditions, but also lifestyle factors, including smoking, hypertension, obesity and excess alcohol consumption (3,4). Individuals who are between 30 and 60 years old experience the highest incidence of IA, and women experience an increased incidence compared with men with a ratio of 3:2 (1,5).

The genes implicated in the pathogenesis of the rupture of IA have been evaluated previously. Through reverse transcription-quantitative polymerase chain reaction, Guo et al (6) identified that the mRNA levels of caspase-3 in IA and abdominal aortic aneurysm are 8.94-fold and 6.73-fold compared with that of normal vessels, which improves the understanding of apoptosis in ruptured intracranial aneurysm. Nuclear factor-κB (NF-κB) functions as an essential regulator during the initiation of IA development via mediation of several inflammatory genes associated with macrophage activation and recruitment, thus, NF-κB may function as a therapeutic target for IA (7). Increased tumor necrosis factor α and Fas-associated death domain protein may have deleterious effects on cerebral arteries through facilitation of inflammation and apoptosis in immune and vascular cells, consequently weakening vessel walls (8,9). Decreased tissue inhibitor of matrix metalloproteinases and increased matrix metalloproteinases in the late stage of IA formation may be the reason for extracellular matrix degradation resulting in the progression, and rupture of IA (10). As a main chemoattractant for monocytes and macrophages, NF-κB activation-induced monocyte chemoattractant protein-1 (MCP-1) acts in IA formation and may serve as a potential target for therapy inhibiting IA progression (11). Nevertheless, the molecular mechanisms of rupture of IA have not been fully revealed.

In 2011, Kurki et al (12) explored the gene expression differences between ruptured and unruptured saccular intracranial aneurysm (sIA) wall samples, screening 686 upregulated and 740 downregulated genes in the ruptured samples, and finding that hypoxia-inducible factor-1A, ETS transcription factors, toll-like receptor signaling and NF-κB are associated with the rupture of sIA walls in humans. In 2010, Pera et al (13) investigated the gene expression profiles of ruptured and unruptured IA, as well as control intracranial arteries, and identified a total of 159 differentially expressed genes (DEGs) and several critical biological processes, including cell adhesion, muscle system, and the immune system, and inflammatory response. However, comprehensive bioinformatic analyses have not been performed to further screen the critical genes associated with the rupture of IA. Using the data gathered by Kurki et al (12) and Pera et al (13), the DEGs between ruptured, and unruptured IA samples were fully screened. Utilizing protein-protein interaction (PPI) network construction, pathway enrichment analysis, construction and efficiency evaluation of support vector machine (SVM) classifier, and principal component analysis (PCA), the key genes implicated in the rupture of IA were further identified.

Materials and methods

Microarray data

Microarray data of GSE13353 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13353) and GSE15629 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE15629) were downloaded from the Gene Expression Omnibus database. GSE13353, which was deposited by Kurki et al (12) and sequenced on the platform of GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array, included 11 ruptured IA wall samples and 8 unruptured IA wall samples. Kurki et al (12) isolated ruptured and unruptured IA wall samples from the necks of Finnish individuals undergoing microsurgical clipping as previously described (14–17). The IA wall samples were frozen in liquid nitrogen and then kept in the Helsinki Neurosurgery sIA Tissue Bank (12). GSE15629 was deposited by Pera et al (13) and sequenced on the platform of GPL6244 (HuGene-1_0-st) Affymetrix Human Gene 1.0 ST Array [transcript (gene) version], from which 8 ruptured IA wall samples and 6 unruptured IA wall samples were selected for this study. Pera et al (13) also resected full-thickness vessel wall samples from patients who underwent microsurgical clipping. Then, the samples were stored in RNAlater (Qiagen China Co., Ltd., Shanghai, China) at 80°C.

Data preprocessing and DEGs screening

Using the oligo package (18) in R, background correction and normalization were performed on the raw data. Using the limma package (http://www.R-project.org) (19) in R, the genes differentially expressed between ruptured and unruptured IA samples were selected, with P<0.05 and |fold change | >1.5 as the thresholds. Conversely, based on the MetaDE package (https://cran.r-project.org/web/packages/MetaDE/index.html) (20) in R, heterogeneity tests and differential expression analysis for each gene were conducted successively. The genes with τ2=0 and Qpval >0.05 were homogeneous and unbiased, from which the genes with P<0.05 were further selected as DEGs.

Construction of PPI network

The interaction information of human proteins was downloaded from the Biological General Repository for Interaction Datasets (http://thebiogrid.org/) (21), Human Protein Reference Database (http://www.hprd.org/) (22) and Database of Interacting Proteins (http://dip.doe-mbi.ucla.edu/) (23) databases. Subsequently, the identified DEGs were mapped to the interaction network of human proteins and the PPI network for the DEGs was visualized using Cytoscape software (version 3.1.0, http://www.cytoscape.org) (24). In the PPI network, nodes and edges separately stand for proteins, and interactions of proteins. Additionally, the number of edges involving one node was the connectivity degree of the node. Furthermore, the candidate genes in the PPI network were further screened using the betweenness centrality (BC) method based on the following formula:

σst represents the number of the shortest path between s and t. σst(ν) stands for the number of paths past node v among the shortest paths between s and t. The value of CB (v) ranges from 0 to 1, a larger value indicates an increased importance.

Additionally, hierarchical clustering analysis (25) was performed for the candidate genes and clustering results were visualized using Heatmap software (version 1.0, http://www.hiv.lanl.gov/content/sequence/HEATMAP/heatmap.html) (26).

Pathway enrichment analysis

Pathway enrichment analysis was performed for the DEGs involved in the PPI network using the fisher algorithm (27) based on the following formula:

Among the formula, M, N and K stand for the number of genes enriched in pathways, the total number of genes in whole genome and the number of DEGs, respectively. p represents the probability that no less than × DEGs are enriched in pathways.

Construction and efficiency evaluation of SVM classifier

As effective classifiers, SVMs may be utilized for two-class classification of microarray data and acquire high classification accuracy (28). The DEGs involved in the PPI network were sorted in descending order based on their BC values. Starting from the last 10 genes and taking 5 genes as an interval, the DEGs were selected as characterization factors (genes used for constructing SVM classifier) in descending order based on BC values. With the dataset of GSE13353 as the training dataset, the tune.svm function of the e1071 package (http://cran.r-project.org/web/packages/e1071/) (29) in R was used to construct an optimal SVM classifier. Then, the dataset of GSE15629 was used as the validation dataset for detecting the SVM classifier. Furthermore, the Candidate Cancer Gene Database (http://ccgd-starrlab.oit.umn.edu/) (30) was used to search for the characterization of gene-associated cancer.

PCA of characterization factors

As a multivariate technique, PCA is able to reassemble associated-indexes into a new set of comprehensive indexes that exhibit no correlations with each other (31). To further optimize the identified characterization genes, the genes involved in the SVM classifier were performed using PCA.

Results

DEGs analysis

In the ruptured samples, a total of 636 DEGs (including 279 upregulated genes and 357 downregulated genes) in GSE13353 and 656 DEGs (including 130 upregulated genes and 526 downregulated genes) in GSE15629 were identified using the limma package. Using the MetaDE package, a total of 1,029 DEGs (including 527 upregulated genes and 502 downregulated genes) were screened. To include the associated genes into our research, the DEGs identified by limma and MetaDE packages were merged and used for the following analyses.

PPI network analysis and pathway enrichment analysis

The PPI network constructed for the merged DEGs identified 510 nodes (including 290 upregulated genes and 220 downregulated genes) and 907 interactions (Fig. 1). Combined with connectivity degrees and BC scores, the top 100 genes [including fibronectin 1 (FN1), amyloid β (A4) precursor protein (APP), nuclear RNA export factor 1 (NXF1) and signal transducer and activator of transcription 3 (STAT3)] in the PPI network were selected as candidate genes. Hierarchical clustering analysis revealed that the candidate genes were able to separate the ruptured samples from the unruptured samples (Fig. 2).

Figure 1.

Protein-protein interaction network for the merged differentially expressed genes. Red and green nodes represent upregulated and downregulated genes in the ruptured samples, respectively.

Figure 2.

Hierarchical clustering trees for (A) GSE13353 and (B) GSE15429. Red and green stand for high and low values, respectively. Rup, ruptured; un, unruptured.

To further understand the biological pathways involved the genes in PPI network, pathway enrichment analysis was conducted. A total of 7 pathways were enriched, including the mitogen-activated protein kinase signaling pathway (P=0.004304), pathways in cancer (P=0.007181) and Toll-like receptor signaling pathway (P=0.040118), which involved toll-like receptor 3 (TLR3) (Table I).

Table I.

Pathways enriched for the genes involved in the protein-protein interaction network.

Description	Gene number	P-value	Gene symbol
hsa04010: MAPK signaling pathway	18	0.004304	FGFR2, ZAK, BRAF, MAP2K3, RELA, MKNK2, MAPK10, STK4, MAX, MAPK1, MAPK12, DUSP1, MAP3K8, SOS2, MAPK8IP2, MAPK8, PAK1, STMN1
hsa05200: Pathways in cancer	20	0.007181	FGFR2, E2F3, RET, BRAF, RELA, STAT5B, ARNT2, EGLN1, MAPK10, STK4, BIRC2, STAT3, MAPK1, MAX, SOS2, RALB, JAK1, CSF3R, MAPK8, FN1
hsa04621: NOD-like receptor signaling pathway	7	0.013053	MAPK1, NOD2, MAPK12, RELA, MAPK8, MAPK10, BIRC2
hsa04012: ErbB signaling pathway	8	0.019615	MAPK1, BRAF, ERBB3, STAT5B, SOS2, MAPK8, MAPK10, PAK1
hsa00270: Cysteine and methionine metabolism	5	0.021434	ADI1, LDHB, LDHA, SRM, TRDMT1
hsa04620: Toll-like receptor signaling pathway	8	0.040118	MAPK1, MAPK12, RELA, MAP2K3, MAP3K8, TLR3, MAPK8, MAPK10
hsa04722: Neurotrophin signaling pathway	9	0.042021	MAPK1, MAPK12, BRAF, RELA, SOS2, NGFRAP1, MAPK8, MAPK10, ARHGDIA

With the dataset of GSE13353 as the training dataset, optimization of the SVM classifier was performed. When the number of characterization genes reached 15, the SVM classifier could completely and accurately distinguish all samples (Fig. 3A). Thus, the 15 genes (including FN1) were considered as the genes required for constructing the SVM classifier (parameters: γ, 0.5; cost, 4; cross, 10). To confirm that the SVM classifier had repeatability and portability, the validation dataset of GSE15629 was used to detect the SVM classifier. In the dataset of GSE15629, the SVM classifier was able to fully differentiate the ruptured samples from the unruptured samples (Fig. 3B). Therefore, the expression pattern characteristics of the 15 characterization genes in ruptured and unruptured samples were notable. Furthermore, the characterization of gene-associated cancer was assessed, revealing that NXF1 was associated with nervous system cancer (Table II).

Figure 3.

Support vector machine classifier was able to fully differentiate the samples in (A) the training dataset of GSE13353 and (B) the validation dataset of GSE15629. Red and black dots represent ruptured and unruptured samples respectively. X and Y axis indicate the position vectors of samples.

Table II.

The 15 characterization genes involved in the support vector machine classifier and the types of cancer associated with them.

Gene	BC_score	Degree	P-value	logFC	Pubmed ID	Cancer type
APP	0.9309	113	0.022006	−1.40543	22057237	Colorectal cancer
BMI1	0.5413	19	0.000794	−1.21654
CCDC8	0.5700	26	0.003424	−1.04234
CFL1	0.5344	14	0.049459	0.817437
FBL	0.5457	17	0.032831	−0.77802
FBXW7	0.5380	19	0.024609	−1.45185	22370638	Blood cancer
FN1	0.6292	42	0.01076	1.299781
HNRNPU	0.5447	24	0.049798	0.906204	27006499	Gastric cancer
NDEL1	0.5339	10	0.036456	0.840461	24316982	Liver cancer
NXF1	0.7250	57	0.004859	1.123627	23685747	Nervous system cancer
RELA	0.5382	18	0.037826	0.860905
STAT3	0.5487	15	0.00697	1.266099	22699621	Pancreatic cancer
TAGLN2	0.5337	12	0.000581	1.101816
TCF4	0.5357	12	0.045941	−0.95848	23045694	Nervous system cancer
UBQLN1	0.5464	16	0.016495	0.942896	24316982	Liver cancer

BC, betweenness centrality; FC, fold change; ID, identification.

To further optimize the identified characterization genes, PCA was conducted on the 15 characterization genes involved in the SVM classifier. APP, NXF1 and STAT3 were the 3 principal components able to separate the samples (Fig. 4).

Figure 4.

Three-dimensional diagram identifying that the three principal components may be able to separate the samples. Red and green circles represent ruptured and unruptured intracranial aneurysm samples, respectively NXF1, nuclear RNA export factor 1; STAT3, signal transducer and activator of transcription 3; APP, amyloid β (A4) precursor protein.

Discussion

In the present study, the DEGs between ruptured and unruptured samples were identified using limma and MetaDE packages. In the ruptured samples, a total of 636 DEGs in GSE13353 and 656 DEGs in GSE15629 were identified by the limma package. Using the MetaDE package, a total of 1,029 DEGs were screened. These DEGs were merged and a PPI network analysis was performed. Combined with connectivity degrees and BC scores, the top 100 genes (including FN1, APP, NXF1 and STAT3) in the PPI network were identified as candidate genes.

Pathway enrichment analysis of the genes in the PPI network revealed that TLR3 was enriched in the TLR signaling pathway. TLRs have been reported to function in vascular inflammatory diseases including aneurysm and atherosclerosis (32). TLR4, which is expressed in IA walls of humans and rats, may be implicated in IA formation via NF-κB activation in endothelial cells (33). FN1 was one of the 15 genes used for constructing the SVM classifier. Alternatively spliced extra domain A (EDA) of fibronectin is essential for tissue repair, and decreased expression of EDA may promote susceptibility to aneurysm of patients with a bicuspid aortic valve (34). Wang and Astrof (35) demonstrated that the local synthesis of FN1 serves essential functions in spatial regulation of Notch signaling and cardiovascular development. As a multi-domain extracellular matrix glycoprotein, fibronectin functions in blood vessel morphogenesis during pathological angiogenesis and embryonic development, and is expressed during pathological angiogenesis in multiple diseases, including late stage artherosclerosis, lung cancer and in abnormal ocular conditions (36–39). Thus, TLR3 and FN1 may serve functions in the rupture of IA.

PCA identified that APP, NXF1 and STAT3 were the three principal components that were able to separate the samples. Duplication of the APP locus, which leads to the increase of amyloid-β peptides, is able to induce autosomal dominant early-onset Alzheimer's disease with cerebral amyloid angiopathy (40). As a result of assessing the characterization of gene-associated cancer, NXF1 was identified to be associated with nervous system cancer. NXF1, which belongs to a family of evolutionarily conserved proteins, includes an NTF2-like domain, a noncanonical RNP-type RNA binding domain, a ubiquitin-associated domain and four leucine-rich repeats (41). In vivo blockade of STAT3 signaling or inhibition of interleukin-17A (IL-17A) lead to an apparent increase in fatal rupture and aneurysm severity in mouse models, and the prevalence of vascular abnormalities are high in patients with STAT3 deficiency (42). Romain et al (43) identified that T cell-specific STAT3 signaling serves a central function in promoting vascular aneurysm and IL-17 serves a protective function in the process. Phosphorylated STAT3 is associated with TLR4-dependent abdominal aortic aneurysm (AAA) formation, and STAT3 and/or TLR4 may be utilized for therapy of AAA (44,45). The level of phosphorylated STAT1 increases during aneurysmal degeneration, and the loss of STAT1 is associated with aneurysm formation and an increased rate of aortic rupture in a model of aortic dissection (46). These findings indicate that APP, NXF1 and STAT3 may be involved in the rupture of IA.

In conclusion, the DEGs between ruptured and unruptured samples were identified using limma, and MetaDE packages. Additionally, TLR3, FN1, APP, NXF1 and STAT3 may function in the rupture of IA. However, a further validation of the roles of these genes in the rupture of IA is required.