Screening of Critical Genes and MicroRNAs in Blood Samples of Patients with Ruptured Intracranial Aneurysms by Bioinformatic Analysis of Gene Expression Data.

BACKGROUND This study aimed to identify more potential genes and miRNAs associated with the pathogenesis of intracranial aneurysms (IAs). MATERIAL AND METHODS The dataset of GSE36791 (accession number) was downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were screened for in the blood samples from patients with ruptured IAs and controls, followed by functional and pathway enrichment analyses. In addition, gene co-expression network was constructed and significant modules were extracted from the network by WGCNA R package. Screening for miRNAs that could regulate DEGs in the modules was performed and an analysis of regulatory relationships was conducted. RESULTS A total of 304 DEGs (167 up-regulated and 137 down-regulated genes) were screened for in blood samples from patients with ruptured IAs compared with those from controls. Functional enrichment analysis showed that the up-regulated genes were mainly associated with immune response and the down-regulated DEGs were mainly concerned with the structure of ribosome and translation. Besides, six functional modules were significantly identified, including four modules enriched by up-regulated genes and two modules enriched by down-regulated genes. Thereinto, the blue, yellow, and turquoise modules of up-regulated genes were all linked with immune response. Additionally, 16 miRNAs were predicted to regulate DEGs in the three modules associated with immune response, such as hsa-miR-1304, hsa-miR-33b, hsa-miR-125b, and hsa-miR-125a-5p. CONCLUSIONS Several genes and miRNAs (such as miR-1304, miR-33b, IRS2 and KCNJ2) may take part in the pathogenesis of IAs.

Background

Intracranial aneurysms (IAs) are a cerebrovascular disease characterized by localized weakened pouches of the arterial wall around bifurcations of the circle of Willis, with disruption of the media and loss of the internal elastic lamina [1,2]. Rupture of IAs can result in subarachnoid hemorrhage (SAH), of which the consequences are dire, carrying an overall mortality rate of about 40% to 50% and a morbidity rate among survivors of approximately 50% [3]. Treatment with both endovascular therapy and microsurgery is effective but can also cause significant morbidity and mortality [4]. Currently, there are no medical alternatives to stabilize aneurysmal progression or prevent rupture. Thus, there is a critical need for a better understanding of the molecular mechanisms underlying the pathogenesis of IAs.

Risk factors predisposing formation of IAs include hypertension, cigarette smoking, heavy alcohol consumption, familial history, cocaine usage, ethnicity, sex, and age [5]. Though the pathophysiology leading to aneurysm formation, growth, and eventual rupture have not been completely elucidated, inflammation is known to play a significant role [6]. Ample evidence has indicated a critical role of TNF-α in the IAs formation and progression to rupture [7,8]. A study demonstrated that hepatocyte growth factor (HGF) concentrations were elevated in blood samples drawn from the lumen of patients with IAs [9]. In addition, Roder et al. performed a meta-analysis of microarray gene expression studies on IAs and revealed that seven genes, including B-cell CLL/lymphoma 2 (BCL2), collagen type I alpha 2 (COL1A2), collagen type III alpha 1 (COL3A1), collagen type V alpha 2 (COL5A2), C-X-C motif chemokine ligand 12 (CXCL12), TIMP metallopeptidase inhibitor 4 (TIMP4), and tenascin C (TNC) were very likely to be involved in the pathogenesis of IAs [10]. On the other hand, the discovery and modulation of microRNAs (miRNAs) provide a novel direction for IA research [11,12]. MiRNAs are small, noncoding RNA molecules of 18 to 25 nucleotides that have fundamental roles in post-transcriptional regulation of gene expression by binding to mRNAs and targeting the mRNA for degradation or translational inhibition [13]. Jiang et al. demonstrated differential expression of 18 miRNAs in IA tissue compared with controls, such as hsa-mir-133b and hsa-mir-143-3p [11]. In addition, the study of Lee et al. showed that miRNAs may take part in cerebral aneurysm formation by affecting multiple target genes and signaling pathways [14]. However, the genes and miRNAs associated with the pathogenesis of IAs were largely unknown.

Recently, Pera et al. performed an analysis of global gene expression profiles in peripheral blood cells from patients with SAH from ruptured IAs and found that rupture of IAs strongly influenced the transcriptional profiles of peripheral blood cells, involving a depression in lymphocyte response with enhancements in monocyte and neutrophil activities [15]. In this study, we downloaded and reanalyzed the dataset of Pera et al. [15] from a publicly accessible database. Blood samples from patients with ruptured IAs and controls were screened for differentially expressed genes (DEGs), followed by functional and pathway enrichment analyses. In addition, gene co-expression network was constructed and significant modules were extracted from the network. MiRNAs that could regulate DEGs in the modules were screened out and the analysis of regulatory relationships was conducted. This study aimed to identify more potential genes and miRNAs associated with the pathogenesis of IAs, which may provide a valuable resource for future biochemical and genetic studies of IAs pathogenesis and rupture.

Material and Methods

Microarray data

The Gene Expression Omnibus (GEO, ) repository is the largest, fully public gene expression resource for high-throughput microarray and other forms of high-throughput data submitted by the scientific community [16]. In this study, we downloaded the microarray data that deposited by Pera et al. [15] from the GEO database [16] with the accession number of GSE36791. The platform is GPL10558, Illumina HumanHT-12 V4.0 expression beadchip. According to the description in the original study of Pera et al. [15], the dataset included peripheral blood samples from 43 patients with SAH from ruptured IAs and 18 controls who were suffering from headaches. SAH was diagnosed using cranial computed tomography and/or lumbar puncture. The presence and location of IAs were evaluated by digital subtractive angiography and/or angio-computed tomography. Venous whole blood was collected before neurosurgical interventions and total RNA was purified from the blood samples. In this study, 43 blood samples from patients with ruptured IAs and 18 blood samples from controls were applied for the following analysis.

Data preprocessing and screening of DEGs

The expression data which had already been normalized using Robust Multi-Array Average (RMA) [17] were downloaded. Then, probe sets were mapped to gene symbols according to the probe annotation files of the GPL10558 platform. When multiple probes mapped to a same gene symbol, the average value of all probes that mapped to the gene was calculated to represent the gene. Afterwards, the gene expression values were log2 transformed.

Linear Models for Microarray Analysis (limma) is an R-based open-source software package that provides an integrated solution for differential expression analyses for data from experiments involving microarrays and other platforms [18]. In this study, the limma package was used to identify DEGs between blood samples from patients with ruptured IAs and controls. False discovery rate (FDR) was estimated using Benjamini-Hochberg method [19]. The DEGs were as defined by FDR <0.05 and fold change (FC) greater than 2 (|log2FC| >0.58).

Gene ontology (GO) and pathway enrichment analysis for DEGs

Database for Annotation, Visualization, and Integrated Discovery (DAVID) bioinformatics resources, is a web-accessible program that enable investigators to discover biological themes associated with large gene lists [20]. In this study, to further characterize the identified DEGs, we performed GO terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the web tool DAVID. Those enriched GO terms and pathways with FDR <0.05 were considered significant.

Construction of gene co-expression network and module extraction

Weighted gene co-expression network analysis (WGCNA) R package provides a comprehensive systems biology method for performing weighted correlation network analysis that describes the correlation patterns among genes and identifying clusters (modules) of highly correlated genes across samples [21]. In this study, in order to identify ruptured IAs-associated co-expression modules and their key constituents, we input expression profiles of the identified DEGs to construct weighted gene co-expression modules using the WGCNA R package [21]. Briefly, the framework for WGCNA involves four steps 1) an unsigned similarity matrix of gene co-expression is initially determined based on pair-wise Pearson correlations for all genes across the samples; 2) the co-expression similarities are transformed into a weighted adjacency matrix of connection strengths using a power adjacency function; 3) a topological overlap matrix that measures the relative interconnectedness (proximity) of the genes is conducted from the adjacency matrix; and 4) average linkage hierarchical clustering of the topological overlap dissimilarity matrix is performed to identify clusters (modules) of co-expressed genes [22,23].

An adjacency matrix (network), A=[aij], was constructed using the following formula:

xi, and xj represent vectors whose components report the gene expression values of gene i and j. cor indicates the Pearson correlation coefficient between xi, and xj. The co-expression similarity matrix was continuously measured with a power (β) as a weight to obtain the weighted adjacency aij. The soft-thresholding parameter b of the power function is selected when the resulting co-expression network (adjacency matrix) best approximates a scale-free topology. In this study, a soft-threshold power of β=6 was used.

Next, the topological overlap-based dissimilarity matrix was computed from the weighted adjacency matrix using the following formula:

where ki indicates the total connectivity of gene i with all other genes in the network.

Average linkage hierarchical clustering was applied using a dissimilarity measure (1-topological overlap matrix) [24,25] and modules which were defined as branches of the resulting hierarchical clustering tree were subsequently identified using the blockwise Modules R function in the WGCNA R package [26]. Each module was subsequently designated by a color. The gene expression profiles of each module were summarized. Moreover, functional enrichment analysis for the DEGs in each enriched module were carried out for identifying overrepresented GO terms and pathways using the DAVID Bioinformatics Resources [20]. A p value <0.05 was set as the criterion.

Screening of disease-related miRNAs and regulatory relationships analysis

The combinatorial Gene Regulatory Networks Builder (cGRNB; ) is a freely available web server to help biologists to model and analyze conditional combinatorial gene regulatory networks through integrated engineering of gene expression datasets and seed-matching sequence information, providing combinatorial regulations involving miRNAs, TFs, and genes [27]. In this study, we downloaded predicted miRNA-gene interactions from cGRNB website [27] with 197,906 regulatory pairs, including 699 mature miRNAs and 8,646 genes. The target genes for a miRNA were defined as Mi. The DEGs enriched in the above modules using WGCNA were defined as N. Fisher’s exact test [28] was applied to measure the significance of the overlap between Mi and N. The algorithm was as follows:

where a is the number of genes that both belong to DEGs in the modules and the target genes of miRNAs; b is the number of genes that belong to the target genes of miRNAs, but not DEGs in the modules; c is the number of genes that belong to DEGs in the modules, but not belong to the target genes of miRNAs; d is the number of genes that neither belong to DEGs in the modules, nor the target genes of miRNAs. In addition, p<0.05 was set as the threshold to identify relevant miRNAs for DEGs in the selected modules. Moreover, the regulatory relationships between the miRNAs and the DEGs in the modules were visualized using Cytoscape software [29].

Results

Data preprocessing and DEGs screening

After data preprocessing, 19,080 genes were obtained. The normalized data for the remaining genes represented similar distribution from sample to sample as shown in the box plot (Figure 1), indicating the reliability of data. Totally, 304 genes were differentially expressed between blood samples from ruptured IAs and samples from controls. There into, 167 DEGs were up-regulated and 137 were down-regulated in blood samples from ruptured IAs.

Figure 1

A box plot of normalized data for differentially expressed genes (DEGs). (A) Data before normalization. (B) Data after normalization.

Functional and pathway enrichment analysis

After enrichment analysis, we found that the up-regulated genes were only enriched in six over-represented GO biological process (BP) terms that were mainly associated with immune response, including inflammatory response, defense response, response to wounding, innate immune response, immune response, and detection of biotic stimulus. By contrast, the down-regulated DEGs were enriched in different GO terms and pathways. Table 1 represented the top 10 GO terms and pathways (ranged by FDR value) that were enriched by the down-regulated genes. The results showed that the down-regulated DEGs were mainly concerned with structure of ribosome and translation.

Table 1

Top 10 gene ontology (GO) terms and pathways enriched by the down-regulated genes.

Category	Term	FDR
GOTERM_BP_FAT	GO: 0006414~translational elongation	1.14E-20
KEGG_PATHWAY	hsa03010: Ribosome	1.11E-15
GOTERM_MF_FAT	GO: 0003735~structural constituent of ribosome	6.80E-14
GOTERM_CC_FAT	GO: 0022626~cytosolic ribosome	2.66E-13
GOTERM_CC_FAT	GO: 0033279~ribosomal subunit	5.33E-13
GOTERM_CC_FAT	GO: 0005840~ribosome	2.24E-11
GOTERM_BP_FAT	GO: 0006412~translation	8.18E-11
GOTERM_CC_FAT	GO: 0044445~cytosolic part	3.14E-09
GOTERM_MF_FAT	GO: 0005198~structural molecule activity	1.64E-07
GOTERM_CC_FAT	GO: 0030529~ribonucleoprotein complex	3.71E-07

The top 10 GO terms and pathways were ranged by FDR value. BP – biological process; CC – cellular component; MF – molecular function; KEGG – Kyoto Encyclopedia of Genes and Genomes; FDR – false discovery rate.

Identification of gene co-expression modules

WGCNA analysis of the DEGs and the following module selection produced six modules significantly altered between ruptured IAs and controls. As shown in Table 2, four out of six modules were enriched by up-regulated genes, whereas genes distributed in two modules were down-regulated. To understand the biological meaning of the modules of up-regulated genes and the modules of down-regulated genes, we carried out functional enrichment analyses. From the results of enrichment analysis, we found that the blue module of up-regulated genes was mainly linked with immune response and anti-apoptosis process. The brown module of up-regulated genes was mainly involved in chromatin assembly and the yellow module of up-regulated genes was correlated with interleukin-1 receptor activity. While the turquoise module of up-regulated genes was concerned with inflammatory response. On the other hand, the blue module of down-regulated genes was mainly enriched in cytolysis and cellular defense response. The turquoise module of down-regulated genes was mainly involved in the translation process. Table 2 only shows the top two GO terms or pathways that were ranked by p values. Collectively, the functions of the blue, yellow, and turquoise modules of up-regulated genes were similar, which were connected with immune response, and had differences with the functions of the brown module of up-regulated genes. The two modules of down-regulated genes were enriched in different functional categories.

Table 2

The top 2 GO terms or pathways enriched by DEGs in the modules.

	Module label	Number of genes	Function	P value
Up_modules	Blue	28	GO: 0000267~cell fraction	0.00278
			GO: 0019864~IgG binding	0.011041
	Brown	21	GO: 0006334~nucleosome assembly	0.0037
			GO: 0031497~chromatin assembly	0.004
	Turquoise	95	GO: 0045087~innate immune response	1.42E-05
			GO: 0006954~inflammatory response	1.00E-04
	Yellow	20	GO: 0004908~interleukin-1 receptor activity	0.0085
			GO: 0019966~interleukin-1 binding	0.012
Down_modules	Blue	16	GO: 0019835~cytolysis	1.70E-04
			GO: 0006968~cellular defense response	0.00159
	Turquoise	119	GO: 0006414~translational elongation	2.30E-25
			hsa03010: Ribosome	5.60E-20

Up_modules represent the modules enriched by up-regulated genes. Down_modules represent the modules enriched by the down-regulated genes.

Screening of ruptured IAs-related miRNAs

Totally, we identified miRNAs for only three modules (blue, yellow, and turquoise) with up-regulated genes of the identified 6 functional modules (Table 3, Figure 2). We found that miR-1304 was significantly enriched with less p value. Among the significant miRNAs, miR-33b could regulate up-regulated genes of two modules (yellow module and turquoise module), such as insulin receptor substrate 2 (IRS2), potassium voltage-gated channel subfamily J member 2 (KCNJ2), and monoamine oxidase A (MAOA). The up-regulated target genes in these two modules for miR-33b ate represented in Figure 3.

Table 3

The predicted miRNAs that regulated the 3 modules of up-regulated genes.

miRNA	Module	P value	Up.DEGs_targets
hsa-miR-1304	Blue	0.003169	CNIH4, CYP1B1, SH3GLB1, YOD1
hsa-miR-373*	Blue	0.004342	PFKFB3, SH3GLB1
hsa-miR-514	Blue	0.012935	ADM, SH3GLB1, YOD1
hsa-miR-33b	Yellow	0.014884	IRS2, KLHL8, MAOA
hsa-miR-568	Turquoise	0.022219	ACSL1, B4GALT5, CREB5, DUSP1, ETS2, GAS7, KLHL2, MEGF9, PHF21A, SIPA1L2
hsa-miR-139-3p	Turquoise	0.023429	BASP1, GAS7
hsa-miR-621	Blue	0.023761	PFKFB3, YOD1
hsa-miR-125b	Turquoise	0.024502	CREB5, FAM126B, KIF1B, MEGF9, MTF1, RNF144B, ZNF281
hsa-miR-575	Yellow	0.033528	IRS2, MANSC1
hsa-miR-125a-5p	Turquoise	0.033995	CREB5, FAM126B, KIF1B, MEGF9, MTF1, RNF144B, ZNF281
hsa-miR-33b	Turquoise	0.03954	FAM126B, KCNJ2, KIF1B, PCNX, RNF144B, ZNF281
hsa-miR-1197	Yellow	0.039541	DAAM2, MAOA
hsa-miR-432*	Turquoise	0.041785	ETS2, RNF144B, ST3GAL4
hsa-miR-1234	Turquoise	0.046918	CEBPB, ZNF281
hsa-miR-545	Blue	0.04737	ADM, CNIH4, SH3GLB1, YOD1
hsa-miR-1208	Yellow	0.048709	DAAM2, IRS2

Figure 2

The miRNAs that regulated the three modules of up-regulated genes. Circular nodes indicate miRNAs. Rectangular nodes denote functional modules.

Figure 3

The up-regulated genes in yellow module and turquoise module targeted by miR-33b. The dotted line denotes the target genes in the yellow module. The solid line denotes the target genes in the turquoise module.

Discussion

In this study, a total of 304 DEGs (167 up-regulated and 137 down-regulated genes) were screened for in blood samples from patients with ruptured IAs compared with those from controls. Functional enrichment analysis showed that the up-regulated genes were mainly associated with immune response and the down-regulated DEGs were mainly concerned with structure of ribosome and translation. Besides, six functional modules were significantly identified, including four modules enriched by up-regulated genes and two modules enriched by down-regulated genes. Therein, the blue, yellow, and turquoise modules of up-regulated genes were all linked with immune response. Additionally, 16 miRNAs were predicted to regulate DEGs in the three modules associated with immune response, such as hsa-miR-1304, hsa-miR-33b, hsa-miR-125b, and hsa-miR-125a-5p.

MiR-1304 is a primate-specific miRNA which has low expression in human brain cortex, peripheral blood, and embryonic stem cells [30,31]. Very few reports are available on the exact roles of miR-1304 to date. Recently, the study of Li et al. demonstrated that miR-1304 suppressed cell growth and survival through inducing cell apoptosis and cell cycle arrest in non-small cell lung cancer [32]. Evidence reports that structural degeneration of aortic tissue with apoptosis of vascular smooth muscle cells (VSMCs) contributes to aneurysmal formation. In the present study, we found that miR-1304 could regulate the up-regulated genes in the blue module with most significance. A single miRNA may regulate several genes as its targets, while one gene may be targeted by many miRNAs. Although the exact role of miR-1304 is unclear, it is reasonable to conclude that miR-1304 may be involved in the pathogenesis of IAs possibly by affecting VSMCs apoptosis.

In this study, miR-33b could regulate the up-regulated genes of two modules (yellow module and turquoise module), such as IRS2, and KCNJ2. Martino et al. indicated that circulating miR-33b were up-regulated in familial hypercholesterolemia in pediatric age cases [33]. In addition, hypercholesterolemia was identified to be risk factor for aneurysms [34]. Thus, we suggested that miR-33b may play a critical role in the development and progression of IAs. In addition, IRS2 encodes the insulin receptor substrate 2, a cytoplasmic signaling molecule which mediates effects of insulin-like growth factor 1, insulin, and other cytokines [35]. A study demonstrated that insulin treatment for glycemic control was safe in patients with aneurysmal SAH [36]. The protein encoded by KCNJ2 is an integral membrane protein and inward-rectifier type potassium channel [37]. Additionally, activation of calcium-activated potassium channels and ATP-sensitive potassium channels may be a major mechanism which regulates vasodilatation of cerebral arteries and arterioles in response to some stimuli [38]. Taken together, these data strongly suggested that miR-33b may be essential in the development and progression of IAs by regulating several genes, such as IRS2 and KCNJ2. However, further studies are warranted to verify these findings.

Our study still had some limitations. First, the sample size was relatively small and further investigations with a larger size of samples are needed. Second, data cross-check and experimental verifications were not included in this study. We will further strengthen our study by further experimental verifications in the future.

Conclusions

This study suggested that several genes and miRNAs (such as miR-1304, miR-33b, IRS2 and KCNJ2) may take part in the pathogenesis of IAs. MiR-1304 may be involved in the pathogenesis of IAs possibly by affecting VSMCs apoptosis. In addition, miR-33b may be essential in the development and progression of IAs by regulating several genes such as IRS2 and KCNJ2. This preliminary work may be a valuable resource, after further characterization in a wider sample set, for designing future studies of both the pathogenesis and pathology of IAs.