This article is mentioned in:

Introduction

Intracranial aneurysm (IA) rupture-induced subarachnoid hemorrhage (SAH) is an acute, non-traumatic and devastating condition that accounts for ~5% of cerebrovascular strokes (1,2), and can result in 45% of mortality, 30% of disability and 50% of cognitive impairment morbidity cases in patients suffering from SAH (3–6). Thus, understanding the molecular mechanism underlying IA rupture-induced systemic damage is critical for public health.

Changes in gene expression profiles in response to disease onset reflect crucial processes involved in disease pathogenesis, development and treatment. The identification of circulating genetic and metabolic biomarkers is of great value for early diagnosis and timely intervention (7,8). For example, a nonsense truncated mutation (R450X, 13:52952757, c.1348 C>T, NM_018676.3) in thrombospondin type 1 domain containing protein 1 (THSD1) is associated with SAH morbidity (9). Furthermore, in zebrafish and mice, THSD1 loss-of-function impaired the focal adhesion of endothelial cells to the basement membrane, resulting in cerebral bleeding and mortality (9). The NF-κB transcription factor (TF), a key factor in inflammation, has been identified to be associated with IA rupture-induced SAH (10,11). In addition, higher expression levels of Toll-like receptor (TLR)4 in the peripheral blood mononuclear cells (PBMCs) of patients with SAH was associated with a larger SAH volume, delayed cerebral infarction and worse functional recovery (12). A therapeutic role for inactivation of TLR4 in neuroinflammation post-SAH has been reported (13). However, these limited reports on recognized genes in response to IA rupture-induced SAH cannot explain the underlying complex neuroinflammation and immune mechanisms post-SAH. IA rupture-induced systemic damage requires improved understanding. The more genes that are changed in response to SAH that can be identified, the clearer the molecular mechanisms underlying the development of cerebral disability and cognitive impairment post-SAH will be.

The present study was conducted using a microarray dataset previously reported by Pera et al (14) with the goal of identifying more circulating biomarkers, including long non-coding RNAs (lncRNAs) and genes in the blood, as well as the pathways associated with them, via systemic and rigorous bioinformatics analyses. Taken together, the study aimed to provide more information on the molecular mechanism underlying SAH-induced cerebral disability and cognitive impairment.

Materials and methods

Data selection

A microarray dataset reporting lncRNA and gene expression profiles in patients with SAH from ruptured IA (GSE36791) (14) was downloaded from the National Center for Biotechnology Information Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). The GSE36791 dataset was generated using the GPL10558 Illumina Human HT-12 V4.0 expression beadchip platform and comprises 4,055 lncRNAs and 19,198 protein-coding genes. The dataset includes 61 samples of venous whole blood, including 43 samples from patients with SAH from ruptured IA (48.0±13.9 years old and 48.8% male) and 18 samples from age-matched healthy controls (46.1±17.5 years old and 72.2% male) suffering from headaches. The study by Pera et al (14) was approved by the ethics committee at their institution and conducted with informed consent from all subjects.

Data processing and gene expression profiling

All txt files in the GSE36791 dataset were downloaded and processed using the Limma package (version 3.34.0; http://bioconductor.org/packages/release/bioc/html/limma.html) (15) for data conversion (log2), background correction and data normalization. The transcript IDs and RefSeq IDs of the lncRNAs and genes in the dataset were annotated in the HUGO Gene Nomenclature Committee database (http://www.genenames.org/) (16). The differentially expressed lncRNAs (DElncRNAs) and genes (DEGs) in the blood samples were screened with the criteria of false discovery rate <0.05 and |log2fold-change (FC)|≥0.263. Hierarchical clustering of the screened DElncRNAs and DEGs in the GSE36791 dataset was performed using pheatmap (version 1.0.8; http://cran.r-project.org/package=pheatmap) (17).

Weighted gene co-expression network analysis (WGCNA) and module identification

To screen lncRNAs and genes with similar expression profiles in patients with SAH from ruptured IA compared with controls, the gene co-expression network and co-expressed gene modules were identified using the WGCNA integrated algorithm (version 1.61; http://CRAN.R-project.org/package=WGCNA) (18,19). All algorithms were conducted following the approximate scale-free features. The co-expression correlation matrix and adjacency function (soft-thresholding parameter and intramodular connectivity) were defined, and the dissimilarity coefficients between genes were calculated. Accordingly, the modules associated with disease traits (phenotypes) were identified and correlation parameters (Pearson's r and P-value) were calculated. Modules of eigengenes (≥50 genes) that are highly associated with disease traits (P<0.05) were identified with the following criteria: cutHeight=0.995 and significant stability correlation (P<0.05). The eigengenes in modules that significantly associated with disease traits were subjected to enrichment analysis.

Enrichment analysis

Gene enrichment was performed to identify the Gene Ontology (20,21) (GO) biological processes (BPs) that were significantly associated with the DEGs in significant modules. The GO BP terms in the Database for Annotation, Visualization and Integrated Discovery (DAVID; v6.8; http://david.ncifcrf.gov/) (22) were picked with the criterion of P<0.05.

lncRNA-mRNA regulatory network construction

The Pearson's correlation coefficients (r) between the DElncRNAs and DEGs in the significant WGCNA modules were calculated using the Cor function (23). The lncRNA-mRNA pairs with r≥0.6 were retained and used as candidates for the construction of a lncRNA-mRNA network. Cytoscape software (version 3.6.1; http://www.cytoscape.org/) (24) was used for the construction of the lncRNA-mRNA network. DEGs in the lncRNA-mRNA network were further exploited for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (25) enrichment analysis with the criterion of P<0.05 for inclusion.

Screening of SAH-associated biomarkers and pathways

Finally, the genes and KEGG pathways associated with IA rupture-induced SAH were screened in the Comparative Toxicogenomics Database (CTD, 2019 update; http://ctd.mdibl.org/) (26) with the key phrase ‘Intracranial Aneurysms’. The overlapping SAH-associated pathways and related genes between the DAVID and CTD databases were selected and regarded as candidate biomarkers or pathways. The DElncRNAs that interacted with the candidate DEGs were selected, and the network involving SAH-associated pathways, DEGs and DElncRNAs was constructed using Cytoscape.

Validation of SAH-associated biomarkers

The expression profiles of the SAH-associated candidates in the GSE36791 dataset were then presented. The expression levels of these DEGs in IA samples and controls were calculated and compared. In addition, the expression profiles of these DEGs in another dataset, GSE13353 (27) [GPL570(HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array], which consisted of 19 wall samples (including 11 ruptured and 8 unruptured IA) were verified.

Statistical analysis

The expression levels of the SAH-associated genes in the GSE36791 and GSE13353 datasets were expressed as the mean ± standard deviation. The differences were identified using a t-test in the Limma package. GraphPad Prism 6.4 software (GraphPad Software, Inc.) was used to generate histograms. P<0.05 was considered to indicate a statistically significant difference.

Results

DElncRNA and DEG statistics in patients with IA

A total of 25 DElncRNAs (13 downregulated and 12 upregulated) and 1,979 DEGs (1,198 downregulated and 781 upregulated) were screened from patients with IA rupture-induced SAH compared with the controls in the GSE36791 dataset (Fig. 1A). The clustering heatmap showed the distinct expression profiles of up- and downregulated DElncRNAs and DEGs in the dataset (Fig. 1B).

Figure 1. Summary of the differentially expressed long non-coding RNA and genes in the peripheral blood samples from patients with intracranial aneurysm (IA). (A) Volcano plot of the DElncRNAs and DEGs. Upregulation and downregulation are presented as blue and red color in the Volcano plot, respectively. (B) Clustering heatmap of the DElncRNAs and DEGs in samples from the GSE36791 dataset. Sample bar: Blue for patients with IA (n=43), and purple for healthy controls with headaches (Con, n=18). Upregulation and downregulation are presented as orange and green color in the clustering heatmap, respectively. FC, fold-change; FDR, false discovery rate; DEG, differentially expressed gene; lncRNA, long non-coding RNA.

Identification of SAH-associated WGCNA modules

WGCNA was conducted according to the approximate scale-free features. Accordingly, an eigengene correlation coefficient square (r2) of 0.9, a soft threshold power of 20, and an intramodular connectivity of 1 were selected as the adjacency function parameters (Fig. 2A). A total of 11 WGCNA modules were identified following the adjacency function parameters in combination with the criteria of eigengenes ≥50 and cutHeight=0.995 (Fig. 2B). Four (purple, turquoise, green and pink) of the 11 WGCNA modules displayed significant positive correlations with disease traits (IA phenotypes; P<0.05 and r>0; Fig. 2C). These modules included 50, 201, 140 and 76 DElncRNAs and/or DEGs, respectively (467 in total; Table SI). One module (black; r=−0.68, n=115) showed a relatively higher negative correlation with disease traits than the gray module (r=−0.67; Fig. 2C). Since most of the DEGs related to IA rupture-induced SAH were upregulated, the present study mainly focused on the selection of upregulated biomarkers; consequently, the downregulated genes in the black module were excluded from analysis. The expression profiles of the upregulated DElncRNAs and DEGs in the four WGCNA modules are shown in Fig. S1.

Figure 2. Adjacency function parameters and construction of WGCNA modules in the GSE36791 dataset. (A) Sample independence and mean connectivity of eigengenes. The red line represents the correlation coefficient square (r2) and mean connectivity of eigengenes under a soft threshold power of 20. (B) Gene dendrogram and WGCNA modules. Gene number >50, cutHeight=0.995, soft threshold power=20 and r2=0.9. (C) Heatmap of the WGCNA module correlations with intracranial aneurysm clinical traits (phenotypes). The blue and red colors indicate negative (r<0) and positive correlations (r>0), respectively. WGCNA, weighted gene co-expression network analysis.

Enrichment analysis for DEGs in SAH-associated modules

In order to identify the biological functions associated with the DEGs in the four SAH-associated modules (purple, turquoise, green and pink), the DEG clusters were individually processed using DAVID, and the DEG-associated GO BP terms were selected. The results indicated that the DEGs in the green, pink and purple modules were associated with GO BP terms related to programmed cell death (including GO:0043067, regulation of programmed cell death; GO:0012501, programmed cell death; GO:0008219, cell death; GO:0012502, induction of programmed cell death; GO:0043068, positive regulation of programmed cell death) and apoptosis (including GO:0043065, positive regulation of apoptosis; GO:0006917, induction of apoptosis; GO:0042981, regulation of apoptosis and GO:0006915, apoptosis). DEGs in the green and purple modules were associated with ‘GO:0009611, response to wounding’, ‘GO:0006952, defense response’ and ‘GO:0010033, response to organic substance’. DEGs in the pink and turquoise modules were associated with ‘GO:0006796, phosphate metabolic process’ and ‘GO:0006468, protein amino acid phosphorylation’. All modules were associated with ‘GO:0006955, immune response’ (Table SII).

lncRNA-mRNA regulatory network construction

All r values between the 467 DElncRNAs and DEGs in SAH-associated modules were calculated, and a total of 1,020 lncRNA-mRNA pairs with r≥0.6 were retained for the construction of the lncRNA-mRNA regulatory network. The network consisted of 1,020 lines (interactions) and 382 nodes (DElncRNAs and DEGs), including 135, 71, 2 and 174 nodes in the green, pink, purple and turquoise module, respectively (Fig. S2). A total of seven DElncRNAs were included, including one green (HCG27, degree=198), one pink [ZNFX1 antisense RNA 1 (ZFAS1), degree=212], no purple and five turquoise [long intergenic non-protein coding RNA (LINC)00265, degree=214; murine retrovirus integration site 1 homolog (MRVI1)-antisense RNA 1 (AS1), degree=128; cytochrome P450 1B1 (CYP1B1)-AS1, degree=115; LINC01347, degree=94; LINC02193, degree=62]. The results demonstrated that mitogen-activated protein kinase (MAPK)-interacting serine/threonine-protein kinase 1 (MKNK1) was regulated by all seven DElncRNAs, whereas Rho GTPase-activating protein 24, BCL6, CXC chemokine receptor 1 (CXCR1), C-terminal domain phosphatase 1, DENN/MADD domain containing 3, NOP2/Sun domain family, member 7 and exportin 6 were regulated by all lncRNAs except LINC02193. To further select the potential hub genes and lncRNAs associated with IA rupture-induced SAH, the lncRNA-mRNA network involving lncRNA-mRNA pairs with r≥0.7 was constructed, which consisted of 197 nodes (190 genes and seven lncRNAs) and 297 lncRNA-mRNA pairs (Fig. 3). Among these pairs, MKNK1 was regulated by four lncRNAs, including MRVI1-AS1, LINC01347, LINC00265 and CYP1B1-AS1. A total of 88 mRNAs, including CXCR1, neutrophil cytosolic factor (NCF)2 and NCF4, were targeted by LINC00265 (Fig. 3).

Figure 3. Regulatory network for the lncRNA-mRNA pairs with high expression correlation. All interactions have a Pearson's correlation coefficient (r) ≥0.7. lncRNAs are represented by boxes and mRNAs are represented by circles. All nodes are upregulated in intracranial aneurysm samples compared with controls. The colors represent the different weighted gene co-expression network analysis modules, and the node size represents the degree of interaction in the network. lncRNA, long non-coding RNA.

KEGG pathway enrichment analysis for DEGs in the lncRNA-mRNA regulatory network

Next, the KEGG pathways associated with all the 375 DEGs in the lncRNA-mRNA regulatory network were screened. DEGs were observed to be enriched in six pathways, including ‘hsa04062: Chemokine signaling pathway’ (CXCR1; CXCR2; VAV3; MAPK1 and α-inhibiting activity polypeptide 3), ‘hsa04060: Cytokine-cytokine receptor interaction’ (CXCR1; CXCR2 and interleukin-4 receptor gene), ‘hsa04010: MAPK signaling pathway’ (MKNK1, MAPK1 and MAPK14), ‘hsa04670: leukocyte transendothelial migration’ (MAPK14, VAV3, NCF2 and NCF4), ‘hsa04620: TLR signaling pathway’ (MAPK1, MAPK14, TLR2, TLR4 and TLR5) and ‘hsa04666: Fc gamma R-mediated phagocytosis’ (MAPK1 and VAV3; Table I).

Table I. KEGG pathways that are associated with differentially expressed genes in the long non-coding RNA-mRNA regulatory network.

Table I.

KEGG pathways that are associated with differentially expressed genes in the long non-coding RNA-mRNA regulatory network.

Term	Count	P-value	Gene names
hsa04062: Chemokine signaling pathway	14	1.35×10−3	VAV3, GNAI3, LYN, PREX1, HCK, RAF1, CXCR1, CXCR2, STAT3, MAPK1, ARRB2, GNB2, PTK2B, GNG5
hsa04060: Cytokine-cytokine receptor interaction	16	4.00×10−3	IL18R1, IL1R2, IL1R1, IL18RAP, CXCR1, TNFSF14, CXCR2, IL17RA, IFNAR1, INHBB, ACVR1B, TNFRSF1A, TNFRSF10C, IL4R, CSF3R, CSF2RA
hsa04010: MAPK signaling pathway	15	1.14×10−2	IL1R2, IL1R1, MKNK1, RAF1, MAPKAPK2, ACVR1B, TNFRSF1A, MAP4K4, MAPK1, ARRB2, MAP3K3, DUSP1, RPS6KA1, RASGRP4, MAPK14
hsa04670: Leukocyte transendothelial migration	9	1.32×10−2	GNAI3, VAV3, NCF2, PTK2B, MAPK14, MMP9, NCF4, VASP, ITGAM
hsa04620: Toll-like receptor signaling pathway	8	1.79×10−2	MAPK1, MYD88, MAPK14, TLR2, FADD, TLR4, TLR5, IFNAR1
hsa04666: Fc γ-R-mediated phagocytosis	7	4.09×10−2	MAPK1, VAV3, LYN, HCK, RAF1, FCGR2A, VASP

Identification of SAH-associated biomarkers and KEGG pathways

The recognized genes and KEGG pathways associated with IA rupture-induced SAH were screened out from CTD using the search phrase ‘Intracranial Aneurysms’. One gene (NCF1) and seven KEGG pathways (Table II) were identified, including three common KEGG pathways between CTD and Table I (hsa04062, hsa04670 and hsa04666; DAVID enrichment analysis). Using the related DElncRNAs, DEGs, and the three common KEGG pathways, an lncRNA-mRNA-pathway network was constructed, which included all the seven DElncRNAs and 21 DEGs [including CXCR2, CXCR1, MAPK1, MAPK14, matrix metalloproteinase 9 (MMP9), NCF2 and NCF4; Fig. 4].

Figure 4. SAH-associated network. SAH-associated KEGG pathways (n=3; white circles), lncRNAs (colored squares) and genes (colored circles) are included. Gray lines indicate interactions between lncRNAs and genes, and red lines indicate links between genes and the KEGG pathways. Potentially important nodes (nodes that were significantly upregulated in intracranial aneurysm with and without rupture and enriched in several KEGG pathways) are presented using larger sizes. SAH, subarachnoid hemorrhage; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNA, long non-coding RNA.

Table II. Intracranial aneurysm (MESH:D002532)-associated genes and KEGG pathways in the Comparative Toxicogenomics Database.

Table II.

Intracranial aneurysm (MESH:D002532)-associated genes and KEGG pathways in the Comparative Toxicogenomics Database.

KEGG pathway	ID	Gene	Overlap
Chemokine signaling pathway	hsa04062	NCF1	Yes
Phagosome	hsa04145	NCF1	No
Osteoclast differentiation	hsa04380	NCF1	No
Fc gamma R-mediated phagocytosis	hsa04666	NCF1	Yes
Leukocyte transendothelial migration	hsa04670	NCF1	Yes
Leishmaniasis	hsa05140	NCF1	No
Fluid shear stress and atherosclerosis	hsa05418	NCF1	No

[i] Overlap indicates the corresponding KEGG pathway is a common pathway between IA (MESH ID, D002532)-associated pathways identified in the Comparative Toxicogenomics Database and the pathways associated with the differentially expressed genes in the blood samples from patients ruptured IA-induced subarachnoid hemorrhage. KEGG, Kyoto Encyclopedia of Genes and Genomes; NCF1, neutrophil cytosolic factor 1; IA, intracranial aneurysm.

Expression profiles of SAH-associated DElncRNAs and DEGs

Finally, the expression profiles of SAH-associated DElncRNAs and DEGs in the GSE36791 dataset were extracted. The GSE13353 dataset, including the gene profiles in ruptured (n=11) and unruptured (n=8) saccular IA wall samples, was also downloaded. The six DElncRNAs (all but LINC01347) and DEGs (including NCF2 and NCF4) in the GSE36791 dataset were all identified in the GSE13353 dataset. The significant upregulation of these eight factors in the ruptured IA blood samples compared with controls in the GSE36791 dataset is shown in Fig. 5A. Significant upregulations of LINC00265, NCF2 and NCF4 were verified in the ruptured IA wall samples compared with unruptured ones in the GSE13353 dataset (Fig. 5B). These findings suggested that LINC00265, NCF2 and NCF4 may serve important roles in post-rupture processes.

Figure 5. Expression profiles of several genes in the GSE36791 and GSE13353 datasets. (A) Expression profiles in the blood samples of controls (n=18) and ruptured intracranial aneurysms (n=43) in the GSE36791 dataset. (B) Expression profiles in the wall samples of patients with ruptured (n=11) and unruptured IAs (n=8) in the GSE13353 dataset. *P<0.05, **P<0.01, ***P<0.001.

Discussion

It has been reported that damage or injury to the cerebral and central nervous system is enlarged and prolonged by SAH-induced inflammation (10,28–30). The inhibition of SAH-induced inflammation attenuates brain injury (31). The present study demonstrated that IA rupture-induced SAH was accompanied by changes in numerous lncRNAs and genes, including HCG2, ZFAS1, LINC00265, MRVI1-AS1, CYP1B1-AS1, LINC01347, LINC02193, CXCR1, CXCR2, TLR2, TLR4, NCF2 and NCF4, as well as pathways, including ‘chemokine signaling pathway’, ‘leukocyte transendothelial migration’ and ‘Fc gamma R-mediated phagocytosis’, which suggested that there are potential roles for these factors in the immune response to SAH.

NCF2 and NCF4 were two upregulated SAH-associated genes. Both NCF2 and NCF4 encode nicotinamide adenine dinucleotide phosphate (NADPH) oxidases (32). Mutations in NCF2, NCF4 and other genes encoding NAPDH oxidases are associated with reduced levels of reactive oxygen species (ROS), higher proportions of neutrophils and incremental perianal disease rate in patients with Crohn's disease (32). Moreover, the NCF2 H389Q mutation predisposes individuals to developing lupus (33), and the mutation in the binding site of NCF4 to phosphatidylinositol 3-phosphate (R58A) promotes ROS reduction and autoimmune arthritis in mice (34). However, NCF2 is induced by tumor necrosis factor-α (TNF-α) (35) and type II interferon-γ (IFN-γ), and is mapped to the IFN signaling pathway (36). The expression of NCF2 is regulated by the TF activator protein 1 (AP-1) by binding to the promoter of NCF2, which could be upregulated by TNF-α (35,37–39). These findings suggest that TNF-α may induce NCF2, which promotes ROS production, activates MMPs (including MMP9 and MMP2) and the p38 MAPK pathway, and enhances inflammation (40–42). The present study identified the upregulation of NCF2 and NCF4, as well as MMP9 and MAPK1, in patients with IA rupture-induced SAH compared with controls and patients with unruptured IA. Thus, the upregulation of NCF2 and NCF4 may indicate enhanced ROS production and inflammatory status in patients with IA rupture-induced SAH.

AP-1 is a switch for various signaling factors (39), which includes NF-κB, TLR, MAPK and JNK (43,44). Activation of ERK, JNK, NF-κB and AP-1, as well as ROS production, in vascular smooth muscular cells could be induced by the pro-inflammatory cytokine TNF-α (44), which has been shown to be elevated in the cerebrospinal fluid of patients with SAH (45). The production of interleukin-6 and TNF-α activates the NF-κB and MAPK pathways, which are associated with an increased percentage of neuronal apoptosis (46–48). It has been reported that suppression of TNF-α, as well as the NF-κB and MAPK pathways, prevents inflammatory arthritis (49), asthma (50) and nervous system injury (46–48). These observations suggest that the upregulation of NCF2, NCF4 and MAPK in patients with IA rupture-induced SAH may indicate their respective roles in SAH-induced inflammation and nervous system injury.

The present study found that NCF2 and NCF4 were regulated by ZFAS1, MRVI1-AS1, LINC00265, HCG27 and CYP1B1-AS1. Among these lncRNAs, LINC00265 has been reported to promote inflammation in abdominal aortic aneurysm (51). LINC00265 knockdown in abdominal artery reduces the formation of aneurysm and inflammation (51). The upregulation of CYP1B1 during inflammation has previously been observed in cancer cells (52), and this link is mediated by the MAPK pathway (53). A study of lncRNA microarray expression profiles in human umbilical vein endothelial cells exposed to high atheroprotective shear stress showed prominent upregulation of CYP1B1-AS1, which was found to have associations with mitotic nuclear division and cell migration (54). lncRNA HCG27 regulates MAPK1 expression by sponging microRNA in gestational diabetes (55). Moreover, ZFAS1 was increased in patients with rheumatoid arthritis, and promoted fibroblast-like synoviocyte migration and invasion (56). However, none of their roles have been reported in inflammation. Their interactions with CXCR1 and CXCR2, two inflammatory chemokine receptors (57–59), indicate roles for these lncRNAs in inflammation processes. Moreover, it has been reported that induction of CXCR1 and CXCR2 results in neuronal death and apoptosis (60,61), whereas the reduction of CXCR1 and CXCR2 attenuates inflammation and neuropathic pain (62,63). These reports suggest that these lncRNAs may serve as therapeutic targets for inflammation post-IA rupture-induced SAH.

In summary, the present study revealed potential roles for lncRNAs and genes in SAH-induced inflammation. The upregulation of lncRNAs (including HCG2, ZFAS1, LINC00265, MRVI1-AS1, CYP1B1-AS1, LINC01347 and LINC02193) and genes (including CXCR1, CXCR2, NCF2 and NCF4) in peripheral blood from patients with IA rupture-induced SAH compared with the controls suggested that they may have potential as biomarkers of SAH or therapeutic targets for inhibiting SAH-induced nervous system injury. However, this hypothesis should be validated using both in vitro and in vivo experiments.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

This study was supported by the National Natural Science Foundation of China (grant no. 81873186).

Availability of data and materials

All data generated or analyzed during this study are included in this published article. The datasets (GSE36791 and GSE13353) analyzed during the current study are available in the the National Center for Biotechnology Information Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/).

Authors' contributions

The study was conceived and designed by XZ and LH. Acquisition, analysis and interpretation of data were conducted by YL, XL and ZC. LH drafted the manuscript. XL revised the manuscript for important intellectual content. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References