Bioinformatics analysis of microRNA profiles and identification of microRNA-mRNA network and biological markers in intracranial aneurysm.

Intracranial aneurysm (IA) is a kind of cerebrovascular disorder, which may result in the subarachnoid hemorrhage with high lethality and disability. The purpose of this study was to reveal the pathogenesis and identify novel biomarkers in IA.We processed the raw microRNA (miRNA) expression profile data of IA obtained from Gene Expression Omnibus. Then weighted correlation network analysis was performed to identify the hub miRNAs in IA. Target genes of hub miRNAs were predicted using multiR package. In addition, a protein-protein network as well as miRNA-mRNA network was constructed and functional and pathway enrichment analyses were done. Finally, the prediction value of hub miRNAs in IA was tested in validation set.Two modules that had relation with IA were identified and 10 hub miRNAs in each module with higher gene-module association were selected. The protein-protein network and miRNA-mRNA network contained 243 nodes and 1496 edges. Functional and pathway enrichment analyses showed that they were mainly enriched in cell cycle, cell proliferation, and PI3K/Akt signaling pathways. Besides, hsa-miR-191-3p, hsa-miR-423-5p, hsa-miR-424-5p, hsa-miR-425-3p were proven to be valuable in prediction IA occurrence.In a word, this study reveals hub miRNAs, target genes and pathways potentially participating in formation and development of IA and screens out some candidate biomarkers. Our findings provide some new perspectives for research and treatment of IA.

Introduction

Subarachnoid hemorrhage is a devastating subtype of stroke affecting younger patients with high mortality and disability. Nearly 85% of cases are due to the rupture of intracranial aneurysm (IA) which results in the increase of intracranial pressure, decrease of cerebral blood flow as well as decrease of cerebral perfusion pressure. It is reported that about 2% to 3% of the general population are affected by IA.[ The diagnosis of unruptured IA relies on the radiological technologies which may cost lots of money and time. However, there is little biomarkers of IA based on peripheral blood test. Furthermore, the specific mechanisms of IA formation and rupture are still largely unknown. MicroRNAs (miRNAs) are the best characterized non-coding RNA, consisting of 18 to 25 nucleotides. These short RNAs are able to inhibit gene expression by binding to the mRNAs. The binding sites of miRNAs are multiple, so that miRNAs can regulate many pathophysiological processes.[ Previous studies have demonstrated that the miRNAs can be detected in circulation including plasma and serum, and these circulating miRNAs can avoid enzymatic degradation via various mechanisms.[ The expression level and type of circulating miRNAs are changeable under different disease conditions suggesting that they play a significant part in the generation and development of diseases and may be served as potent diagnosis biomarkers. Lots of studies have proven that miRNAs are important in modulating the cell proliferation, apoptosis, inflammation so that they can influence many diseases such as myocardial infarction, hypertension, cancer even aneurysms.[ However, we know little about the function of circulating miRNAs in aneurysms especially in IA. Compared with the miRNAs localized in tissue, circulating miRNAs are more convenient for detecting and have a wider range of influence.

The invention of novel algorithm including weighted correlation network analysis (WGCNA),[ gene set enrichment analysis[ and establishment of database such as gene ontology (GO) knowledgebase,[ Kyoto encyclopedia of genes and genomes (KEGG) database[ make it possible to explore the pathophysiology mechanisms of diseases in a more comprehensive and profounder way. Therefore, in this study, we will use the bioinformatics methods combined with public databases to comprehensively explore the potent function of circulating miRNAs in IA and screen out the hub miRNAs that may be considered as potential biomarkers.

Materials and methods

Microarray data processing

The raw data of a human IA miRNAs microarray dataset (GSE50867) were downloaded from GEO.[ The R package Limma was used to process the raw data.[ Probe IDs were converted into miRNA names according to the proper GPL file downloaded from GEO. To facilitate the further analyses, we transformed the miRNA names into latest version according to the miRBase (v22.0).[ Expression clusters and box plots were done to check whether there was a batch effect or v outliers.

Construction of co-expression networks

The WGCNA R package was used to contrast the co-expressed gene networks. Modules were generated using the dynamic tree-cutting function with following parameters: TOMType = “signed,” deepSplit = 2, pamStage = TRUE, cutHeight = 0.99, and minClusterSize = 30-3∗deepSplit. (The module preservation function (nPermutations = 200) of the WGCNA package (Langfelder et al, 2011) was utilized. Modules showing 75% similarity are merged together. Module eigengene were correlated with clinical traits and the correlations were calculated by linear regression model. The higher correlation means the more closely connection between modules and clinical traits. Hub miRNAs are those that show higher gene-module association KME values.

Gene set variation analysis (GSVA) score

To verify the high correlation between hub miRNAs and disease state, we used R package GSVA to calculate the score of each control sample or IA sample.[

Target mRNA identification and network construction

Since miRNAs play a biological role by modulating the mRNA, we firstly identified target mRNA of hub miRNAs using R package multiMiR.[ This package consist of multi miRNAs databases and all data can be divided into 2 parts according that the targets of miRNAs are validated or predicted. The validated data can be obtained from mirecords, mirtarbase and tarbase, and the predicted data can be acquired from diana_microt, elmmo, microcosm, Miranda, mirdb, pictar, pita as well as targetscan. The string database was used to explore the protein-protein interaction (PPI) of target, and cytoscape was used to construct and show the miRNAs-mRNAs network.[

GO and KEGG pathway analysis of gene expression networks

The validated and predicted target gene gained via multiR were processed utilizing the R package clusterprofile.[ We tried to annotate they biological function via GO and KEGG databases. The output combined scores and adjusted P-values were used to determine significance.

Validating the prediction performance of hub miRNAs

To verify the prediction performance of hub miRNAs, another independent expression set downloaded from GEO (GSE66239) was used as validation set.[ Receiver operating characteristic (ROC) curves were constructed using R package pROC.[

Results

Raw data processing and quality control

The raw data of GSE50867 were downloaded from GEO and handled with the Limma package. The box plot of processed data suggested that there is no obvious bias among each group, and the expression clustering results showed no outlier samples. Principal component analysis shows significant segregation of different group samples. All these results are shown in Figure 1.

Figure 1

Quality control of processed data. (A) Boxplot of processed data. (B) Expression clusters of processed data. (C) Principal component analysis of processed data.

Weighted correlation network analysis (WGCNA) of processed data. (A, B) Scale independence and mean connectivity analysis. (C) Cluster dendrogram among modules. (D) Module-trait relationships. The darker the module color, the more significant their relationship.

WGCNA and hub gene identification

To analyze relationships between miRNAs and IA, the entire data were employed to constructed a weighted gene co-expression network which clusters miRNAs based on their expression in all investigated samples into co-expression modules. According to the Figure 1A, a scale-free network distribution with higher connectivity was formed when the soft-threshold power β was set to 5. As we can learn from Figure 1C, brown module showed the highest positive correlation with IA (the P-value was .0088), meanwhile, the green module had the most significant negative correlation (the P-value was .0039). Thus, we further screened the top 10 hub miRNAs in 2 modules based on the module membership (kME) values. hsa-miR-363-3p, hsa-miR-192-5p, hsa-miR-425-5p, hsa-miR-25-3p, hsa-miR-423-5p, hsa-miR-20b-5p, hsa-miR-103a-3p, hsa-miR-424-5p, hsa-miR-140-5p, hsa-miR-93-5p are 10 hub miRNAs in brown module and hsa-miR-1281, hsa-miR-1825, hsa-miR-498-5p, hsa-miR-1280, hsa-miR-1234-3p, hsa-miR-720, hsa-miR-1228-3p, hsa-miR-425-3p, hsa-miR-191-3p, hsa-miR-1225-3p are 10 hub miRNAs in green module. All these miRNAs may have major impact on the generation and development of IA, and may regarded as the potential biomarkers. To verify the reliability of these hub miRNAs, we referred them to 2 different gene sets according their affiliation and used GSVA to calculate the scores. As shown in Figure 3A, the score of brown module in control was significant lower than that in other 2 groups (scores: control group, −0.78 ± 0.10; unruptured group, 0.65 ± 0.11; ruptured group, 0.44 ± 0.11, P < .05, 1 way ANNOVA). Meanwhile, Figure 3B showed that compared with control group, the scores of green in IA groups were decreased. (scores: control group, 0.75 ± 0.05; unruptured group, −0.55 ± 0.17; ruptured group, −0.52 ± 0.13, P < .001, 1 way ANNOVA). A conclusion drawn from these 2 distinct algorithms was that hub miRNAs identified above were highly related to IA and may play an extremely important role in the pathophysiology of this disease. What's more, considering their great correlation with IA, they may be served as forecast indicators.

Figure 3

Gene set variation analysis (GSVA) scores. (A) GSVA score of brown module hub miRNAs in GSE50867. (B) GSVA score of green module hub miRNAs in GSE50867. miRNAs = microRNAs.

Target mRNA identification and enrichment analysis

Since miRNAs work through regulating mRNAs, to further understand the function of hub miRNAs in IA, the multiMiR was employed to identify the target mRNAs. MultiMiR package provides 2 kinds of databases including validated and predicted data. To filter out the mRNAs that may have huge influence on IA, we included targets which are modulated by more than 3 disparate miRNAs, then enrichment analysis was done to reveal their function. When it comes to the target genes of green module hub miRNAs, there is no terms can be enriched through GO or KEGG analysis, consequently, we paid attention to the brown module hub miRNAs in the following. As shown in Figure 4A, when it referred to genes controlled by brown module hub miRNAs, GO biological process terms were mainly associated with tissue development, cell cycle, cell proliferation, cell division, cell adhesion. In the molecular function group, GO terms mainly involved promoter binding, transcription activator or coactivator binding (Fig. 4B). There are no items enriched in cellular components (CC) and KEGG pathway.

Figure 4

Gene Ontology (GO) enrichment analysis of brown module hub miRNAs target mRNAs. (A) Biological process (BP). (B) Molecular function (MF). miRNAs = microRNAs.

PPI and miRNA-mRNA network

With the help of STRING online database (https://string-db.org), a PPI network was constructed to study the relationship among each mRNA and cytoscape was utilized to build and visualized the miRNA-mRNA network.[ The network containing 233 different genes and 10 hub miRNAs, connected with 1496 edges (Fig. 5 A). Among the 243 nodes, the most significant 10 node degree genes were PTEN, VEGFA, CCND1, MDM2, CREB1, AGO2, BCL2L11, TNRC6B, SMAD7, and KIF23, and the hub miRNA that connected with others most extensively is hsa-miR-93-5p. To further understand the purposes of this network, annotation was performed. As we can conclude from Figure 5 B, PI3K-Akt signaling pathway, FoxO signaling pathway, p53 signaling pathway as well as other cancer related pathways were significantly enriched, suggesting the involvement of cell proliferation, cell adhesion. Besides, another important pathway named (mitogen-activated protein kinase) MAPK was also enriched, indicating that miRNAs may extensively participate in IA generation and development. CC analysis showed that this PPI network is strongly correlated with cell adhesion, different junction, and extrinsic component of membrane (Fig. 5 E).

Figure 5

Protein-protein interaction (PPI) and miRNA-mRNA network and enrichment analysis. (A) The bigger of label size, the more significant of node in network. The hexagon represents the miRNAs and the circle represents the target mRNAs. (B) Signaling pathways enrichment of PPI network. (C) BP. (D) MF. (E) Cellular component (CC). BP = biological process, MF = molecular function, miRNAs = microRNAs.

Expression level and prediction performance of hub miRNAs

To validate the prediction performance of hub miRNAs, another independent dataset was utilized. All identified hub miRNAs belonging to green or brown module were used respectively to draw the ROC curve with the assistance of R package (pROC). We chose the hub miRNAs whose area under curve (AUC) value is higher than 0.75. As shown in Figure 6, the hsa-miR-191-3p performs better than others with the AUC being 100%. Meanwhile, hsa-miR-423-5p (AUC 75.9%), hsa-miR-424-5p (AUC 83.3%), hsa-miR-425-3p (AUC 75.9%) also had good prediction performance. In order to clarify their potential biomarker roles more clearly, we show the expression level of hub miRNAs in Figure 6E.

Figure 5 (Continued)

Protein-protein interaction (PPI) and miRNA-mRNA network and enrichment analysis. (A) The bigger of label size, the more significant of node in network. The hexagon represents the miRNAs and the circle represents the target mRNAs. (B) Signaling pathways enrichment of PPI network. (C) BP. (D) MF. (E) Cellular component (CC). BP = biological process, MF = molecular function, miRNAs = microRNAs.

Receiver operating characteristic (ROC) curves for hsa-miR-423-5p (AUC 75.9%), hsa-miR-424-5p (AUC 83.3%), hsa-miR-425-3p (AUC 75.9%), hsa-miR-191-3p (AUC 100.0%). AUC indicates area under the curve. (E) The expression level of hub miRNAs. miRNAs = microRNAs.

Discussion

IA is a kind of cerebrovascular disorder which mainly performs the abnormal dilatation or ballooning in the blood vessel, localized in the thickness part of the vascular wall. The specific mechanism of its generation and development is still unclear and little potent biomarkers are identified. Recently, circulating miRNAs are viewed as potentially valuable predictors for their stability, sensitivity as well as broadly influence.[ Furthermore, previous studies have demonstrated that miRNAs are significant in regulating the IA.[ Thus, in this studies, with the assistance of WGCNA, we obtained the modules highly associated with IA and screened twenty key miRNAs. Then, the GO terms and pathway enrichment analysis of target mRNAs or PPI network were enriched in biology associated with cellular characteristics, especially for the cell cycle, cell proliferation, cell adhesion, and promoter binding. Many former investigations have demonstrated that the apoptosis and proliferation of endothelial cell and smooth muscle cell, immune cells infiltration are 2 major causes of IA.[ Our results further prove that these 2 main pathophysiological processes can be adjusted by miRNAs and indicates that miRNAs play an important role in IA.

According to our results, PI3K/Akt signaling pathway is important in IA which is consistent with previous studies. Li et al have proven that PI3K/Akt signaling pathway has effect on the production of IA in rat model.[ The involvement of this pathway in other kinds of aneurysm such as abdominal aortic aneurysm and thoracic aortic aneurysm is also detected by different groups.[ Based on existing knowledge, we speculate that PI3K/Akt may affect IA via following ways. First, it may regulate the tight junction of endothelial cells, resulting in the immune cells infiltration and intimal structure changes.[ Second PI3K/Akt may modulate the proliferation and migration of vascular smooth muscle cell (VSMC).[ The balance between apoptosis and proliferation of VSMCs as well as the correction arrangement are crucial in IA growth. Thirdly, PI3K/Akt pathway can control the activation of NF-κB via IκB, published data support that NF-κB essential in regulating the inflammation.[ Dysregulated inflammation exacerbates the IA via various ways, for example, degrading the extracellular matrix through matrix metalloproteinase and recruiting the immune cells.[ What's more, FOXO1, one of the downstream transcription factor of PI3K/Akt pathway, has been proven to be capable of binding with the promoter region of CCND1. This process can modulate the proliferation of VSMCs and vascular remodeling, which influence the development of IA.[

Besides the PTEN which is vital molecule in above mentioned PI3K/Akt pathway, VEGFA also takes up a core position in the PPI network. Li et al have verified that when binding with its receptors, VEGFA contribute to abdominal aortic aneurysm formation by multiple ways, for instance, suppressing mural angiogenesis, producing matrix metalloproteinase and inducing myeloid cell as well as monocytes chemotaxis.[ The similar functions of VEGFA in IA have also been found in IA by a different group.[

Serving circulating miRNAs as biomarkers have attracting considerable interest and pioneering studies have revealed the changes in plasma miRNAs in response to IA.[ In this study, we used WGCNA and identified 4 miRNAs which might be potential biomarkers in IA, with AUC value higher than 75%. In addition to concerning the sensitivity and specificity, the expression levels of latent biomarkers are also crucial, because of the possibility of miRNAs degradation. All these 4 miRNAs are abundant in samples. Consequently, our results suggest that they may be biological markers that are useful in evaluation the probability of IA occurrence.

However, there are still many limitations in our research. First of all, we only used the data obtained from GEO, but we do not check this result in our specimens. Second, the total number of sample is 12 (4 per group), which might be too little. Third, the functions of hub miRNAs predicted through databases are needed to be further verified using cell experimentation and clinical samples. At last, the validation set and training set come from different organizations. Therefore, the results, especially the predictive performance of hub miRNAs could be influence by samples derived from different sources.

In conclusion, we identified the key miRNAs, their target mRNAs as well as crucial pathways. What's more the possibility of using hub miRNAs as biomarkers are verified. These findings are beneficial to comprehensive understanding of potential mechanism in IA and provide some novel perspectives on the IA therapeutic targets.

Author contributions

Data curation: Ming Zhao, Longbiao Xu.

Formal analysis: Ming Zhao.

Funding acquisition: Ming Zhao, Hui Qian.

Investigation: Hui Qian.

Project administration: Ming Zhao.

Software: Ming Zhao, Longbiao Xu.

Visualization: Ming Zhao, Longbiao Xu.

Writing – original draft: Ming Zhao, Longbiao Xu.

Writing – review and editing: Ming Zhao, Hui Qian.