Skip to main content

Comprehensive analysis of ceRNA network related to lincRNA in glioblastoma and prediction of clinical prognosis

Abstract

Background

Long intergenic non-coding RNAs (lincRNAs) are capable of regulating several tumours, while competitive endogenous RNA (ceRNA) networks are of great significance in revealing the biological mechanism of tumours. Here, we aimed to study the ceRNA network of lincRNA in glioblastoma (GBM).

Methods

We obtained GBM and normal brain tissue samples from TCGA, GTEx, and GEO databases, and performed weighted gene co-expression network analysis and differential expression analysis on all lincRNA and mRNA data. Subsequently, we predicted the interaction between lincRNAs, miRNAs, and target mRNAs. Univariate and multivariate Cox regression analyses were performed on the mRNAs using CGGA data, and a Cox proportional hazards regression model was constructed. The ceRNA network was further screened by the DEmiRNA and mRNA of Cox model.

Results

A prognostic prediction model was constructed for patients with GBM. We assembled a ceRNA network consisting of 18 lincRNAs, 6 miRNAs, and 8 mRNAs. Gene Set Enrichment Analysis was carried out on four lincRNAs with obvious differential expressions and relatively few studies in GBM.

Conclusion

We identified four lincRNAs that have research value for GBM and obtained the ceRNA network. Our research is expected to facilitate in-depth understanding and study of the molecular mechanism of GBM, and provide new insights into targeted therapy and prognosis of the tumour.

Peer Review reports

Background

Glioblastoma (GBM) is the most aggressive and destructive primary malignant central nervous system tumour [1]. At present, the most established treatment protocol for this tumour involves curative surgery and adjuvant radiotherapy combined with temozolomide; however, its prognosis is still very poor, and the average survival time of patients is approximately 2 years [2]. Moreover, various treatment methods affect normal brain tissue, and the quality of life of patients does not significantly improve afterwards. Therefore, an in-depth study of the molecular mechanism of GBM occurrence and development, while exploring possible therapeutic targets, may prove beneficial for the diagnosis and treatment of the tumour.

Although approximately 90% of genes in the human genome can be transcribed, only approximately 2% of the genes are protein coding, and non-coding RNAs account for most of the remainder [3]. Long non-coding RNAs (lncRNAs) contain > 200 nucleotides and have no protein coding function. They are mainly transcribed from different regions of the genome by RNA polymerase II, and have been shown to be closely related to cancer [4, 5]. The tissue specificity of lncRNA expression has been reported to be higher than that of mRNA expression, which applies even in pathological conditions such as cancer [6]. Long intergenic non-coding RNA (lincRNA) is a type of lncRNA that does not overlap with exons of protein-coding genes and other non-lincRNA genes, and participates in many important biological processes [7, 8]. Cancer-related lincRNAs may be targeted to provide new ways for cancer diagnosis and treatment [9]. Linc-ROR is a tumour promoter, which mainly participates in tumour cell proliferation, apoptosis, invasion, angiogenesis, and cancer stem cell generation by regulating target genes [10]. NEAT1 is a p53-regulated lincRNA, which plays a key role in cancer occurrence [11]. LincRNA-p21 plays an important role in regulating TAM function in tumour microenvironments, and can be used as a new cancer treatment target for macrophage infiltration [12]. In general, molecular mechanisms related to lincRNA in tumours have important research value.

The COX proportional hazard regression model can analyze the impact of several factors on survival simultaneously, and is the most commonly used multi-factor analysis method for survival analysis. Long et al. found that a prognostic model composed of 4 mRNAs is a reliable tool for predicting the survival of liver cancer patients [13]. A new type of COX model containing 6 genes can improve the prognosis prediction of patients with uterine sarcoma [14]. In summary, the COX model can help clinicians select personalized treatment plans for patients and predict prognosis.

A competitive endogenous RNA (ceRNA) network is the interaction between lncRNAs, miRNAs, and mRNAs that constitutes a complex regulatory network system, which has extensive functions in the human genome and plays a significant role in cancer [15]. The ceRNA network includes mRNAs, miRNAs, lncRNAs, and circRNAs, which play a key role in the occurrence and development of gastric cancer and colorectal cancer, and can be used as biomarkers for tumour treatment [16, 17]. It has been reported that linc-ROR can be used as the ceRNA of miR-145 to regulate the proliferation, invasion, and tumourigenicity of pancreatic cancer cells [18]. Further, a study has demonstrated that lincRNA-HOTAIR is highly expressed in a variety of tumour tissues and cells, and is associated with tumour metastasis and poor prognosis. Additionally, its related ceRNA network has been shown to play a role in the progression of kidney cancer [19]. At present, there are few studies on the ceRNA network in glioblastoma, conducted mainly for a single database and with a small sample, and there are fewer prognostic prediction analyses [20, 21]. In their study, Liu et al. used two databases to construct a ceRNA network for GBM, but did not incorporate WGCNA analysis, nor did they systematically construct a prognostic prediction model [22]. Therefore, in this study, we comprehensively analysed the GBM data retrieved from The Cancer Genome Atlas (TCGA), National Center for Biotechnology Information (NCBI), Gene Expression Omnibus (GEO), Chinese Glioma Genome Atlas (CGGA), and Genotype Tissue Expression (GTEx) databases. We used weighted gene co-expression network analysis (WGCNA) and differential analysis to screen key lincRNAs, construct lincRNA-related ceRNA networks in GBM, and obtain a prognostic model. We believe this study provides a basis for further studies on the molecular mechanism of GBM and exploration of new therapeutic targets.

Methods

Data acquisition and pre-processing

Figure 1 shows the flow diagram of our study. We obtained RNA-Seq data of 165 GBM samples and 1029 normal brain tissue samples from the TCGA and GTEx databases, respectively, using the UCSC Xena Browser (https://xena.ucsc.edu/). Further, the mRNA-seq and clinical data of 693 glioma samples were retrieved from CGGA (http://cgga.org.cn/), from which 248 GBM samples were selected for this study. Among the selected samples, 236 patients had survival information, while 198 had all clinical information. Additionally, we retrieved GSE50161 (comprising 34 GBM and 13 normal samples; platform: Affymetrix-GPL570), GSE134783 (comprising 71 GBM samples; platform: Affymetrix-GPL570), and GSE90603 (comprising 16 GBM and 9 normal samples; platform: Affymetrix-GPL21572) from the GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Fig. 1
figure1

Study workflow. TCGA, The Cancer Genome Atlas; GTEx, Genotype Tissue Expression; GEO, Gene Expression Omnibus; CGGA, Chinese Glioma Genome Atlas; GBM, glioblastoma; WGCNA, Weighted Gene Co-expression Network Analysis; DEG, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis

The IDs of all samples were transformed into gene symbols based on GENCODE (https://www.gencodegenes.org/human/). The RNA-seq data fetched from TCGA and GTEx were collated and merged, and the expression matrices of mRNA and lincRNA were obtained. We also combined the GEO data and obtained the mRNA and miRNA expression matrix. The unit of RNA-seq data is FPKM. The R package oligo [23] was utilized for format conversion, missing data filling, background correction, and data normalisation. All data were downloaded from public databases, and the data were approved by medical ethics when they were first published. Therefore, this study did not require medical ethics documents.

Construction of weighted gene co-expression network

Using the combined data of TCGA and GTEx, a weighted gene co-expression network was constructed for lincRNA and mRNA. Based on the mRNA of the GEO combined data, a weighted gene co-expression network was constructed. WGCNA R software package [24] was used to build a co-expression network to mine key modules related to GBM. First, a hierarchical cluster analysis was performed to check the heterogeneity of the samples. Subsequently, according to the scale-free topological standard, we chose the appropriate soft threshold power (β) to construct a weighted adjacency matrix and convert the adjacency relationship into a topological overlap matrix. Finally, we obtained the gene modules, and determined the modules related to the pathogenicity of GBM through the relationship between the modules and traits, and used a scatter plot to show the correlation between gene significance and module membership.

Screening of differentially expressed genes

We used the Limma R package [25] to screen the differentially expressed lincRNA (DElincRNA) between the GBM patients and control group from the combined data of TCGA and GTEx.

Differential expression analysis was performed on the mRNA data in the TCGA and GTEx combined data and GEO combined data, and the two selected differential genes were combined to obtain differentially expressed mRNA (DEmRNA). Additionally, differential analysis was performed on GSE90603 to obtain DEmiRNA. (P < 0.05, and |log2 FC| ≥ 1).

Preliminary construction of ceRNA network

The most relevant modules to GBM were selected and analysed with DElincRNA to screen for overlapping lincRNA. Further, miRcode (http://www.mircode.org/) was used to predict miRNA, and the miRNA target mRNA was predicted using miRDB (http://www.mirdb.org/), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), and TargetScan (http://www.targetscan.org/). The predicted mRNA, WGCNA key module mRNA (TCGA and GTEx, GEO), and DEmRNA were analysed to obtain the overlapping mRNA. Additionally, Cytoscape 3.7.2 software was used to visualise the ceRNA network.

Gene function and pathway enrichment analysis

We used the clusterProfiler R package [26] for Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. GO is used to describe gene function, and KEGG is used to obtain possible pathways. Further, the GOplot R package [27] was used to visualise the GO term or KEGG approach (adjusted p < 0.05).

Construction of cox regression model

Univariate COX regression analysis was performed on the CGGA data, using the Survival R software package, to evaluate the effect of mRNA expression on the survival time of GBM patients. Additionally, we performed multivariate COX regression analysis, constructed a multivariate Cox regression model, and identified the corresponding coefficients of GBM prognostic features. Further, we calculated the risk score to predict survival time, dividing the samples into high-risk and low-risk groups, with the median as the critical value. The “predict ()” function was used to calculate the risk score: risk score = h0(t)*exp.(β1x1 + β2x2+ … + βnxn). The correlation between the prognostic characteristics of patients and overall survival rate was then calculated through univariate and multivariate Cox regression analyses of clinical factors related to overall survival. Finally, the R package survivalROC was used to draw the receiver operating characteristic curve and calculate the area under the curve (AUC).

CeRNA network in GBM

Based on the mRNA obtained through multivariate Cox regression analysis and DEmiRNA, we screened and reconstructed the ceRNA network.

Gene set enrichment analysis (GSEA)

GSEA software was used to analyse the RNA-seq data of lincRNA retrieved from the CGGA database. According to the median of lincRNA expression, the samples were divided into high and low expression groups. Statistical significance was set at p < 0.05. A set of genes “c2. cp. kegg. v7.2. symbols. gmt” was retrieved from the Molecular Signature Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp) and selected as the reference gene set.

Results

WGCNA identifies key modules

We combined GTEx and TCGA data to construct a co-expression network for all lincRNAs, using the R package ‘WGCNA’, and confirmed that the β value in the network was 4 (Fig. 2a). Further, the dynamic tree cutting method was used to generate co-expression modules, and the closely related modules were merged into larger modules. Finally, 19 modules were generated in the lincRNA co-expression network (Fig. 2b). The module eigengenes (MEs) of the turquoise module had the strongest correlation with GBM traits (Fig. 2c). Figure 2d shows the correlation between gene significance and module membership in the turquoise module, which was considered a key module containing 2206 lincRNAs.

Fig. 2
figure2

Identification of lincRNAs modules highly related to traits. a, Analysis of network topology for various soft-thresholding powers. Left: Analysis of the scale-free fit index for various soft-thresholding powers (β). Right: Analysis of the mean connectivity for various soft-thresholding powers. b, Cluster dendrogram of all lincRNA in the co-expression network. As a result, 19 co-expression modules were constructed and was shown in different color. c, Module-trait associations of lincRNAs were evaluated by correlations between MEs and GBM pathogenicity. Each row corresponds to a module eigengene, column to a trait. Each cell contains the corresponding correlation and p-value. d, The correlation between Gene Significance (GS) and Module Membership (MM) in the turquoise module. There is a highly significant correlation between GS and MM

We also constructed the co-expression networks of all mRNAs for GTEx and TCGA combined data and GEO combined data, comprising 20,270 and 14,807 mRNAs, respectively. In the mRNA co-expression networks of the two sets of data, β values were 8 and 4 (Fig. 3a and Fig. 4a), and 22 and 21 modules were generated (Fig. 3b and Fig. 4b), respectively. In the combined GTEx and TCGA data, the MEs of the blue and dark green modules have the strongest correlations with the tumour and normal traits, respectively (Fig. 3c). The blue and dark green modules were considered key modules, comprising 7736 mRNAs. In the GEO data, Fig. 4c shows that the MEs of the blue and turquoise modules (comprising 8022 mRNAs) have the strongest correlations with the tumour and normal traits, respectively. Figures 3d-e and Fig. 4d-e show the correlation between gene significance and module membership.

Fig. 3
figure3

Identification of mRNAs modules highly related to traits in the combined data of TCGA and GTEx. a, Determination of soft-thresholding power. b, Cluster dendrogram of all mRNA in the co-expression network. c, The correlation between modules and GBM pathogenicity were displayed. d, The correlation between GS and MM in the blue and darkgreen modules

Fig. 4
figure4

Identification of mRNAs modules highly related to traits in the combined data of GEO. a, Determination of soft-thresholding power. b, Cluster dendrogram of all mRNA in the co-expression network. c, The correlation between modules and GBM pathogenicity were displayed. d, The correlation between GS and MM in the blue and turquoise modules

Identification of differentially expressed lincRNAs (DElincRNAs), miRNAs (DEmiRNAs) and mRNAs (DEmRNAs)

We identified DElincRNAs and DEmRNAs from the GTEx and TCGA combined data, including 163 up-regulated lincRNAs, 176 down-regulated lincRNAs, 2953 up-regulated mRNAs, and 2932 down-regulated mRNAs. From the GEO data, 2434 up-regulated and 1995 down-regulated mRNAs were identified. The two groups of mRNAs were then comprehensively analysed to obtain overlapping 1040 down-regulated mRNAs and 1358 up-regulated mRNAs. Finally, 339 DElincRNAs and 2398 DEmRNAs were identified. Through the GSE90603 dataset, 155 down-regulated and 134 up-regulated miRNAs were identified.

Preliminary construction of ceRNA network

Through WGCNA analysis of lincRNAs, we obtained 2206 lincRNAs in the turquoise module and then integrated the analysis with 339 DElincRNAs to obtain 251 overlapping lincRNAs. Further, using the miRcode online tool, 251 lincRNAs were used to predict miRNAs. Additionally, using the miRDB, TargetScan, and miRTarBase datasets, the predicted miRNAs were used to obtain the corresponding mRNAs. By analysing the predicted mRNAs, WGCNA key module mRNA (TCGA and GTEx, GEO), and DEmRNA, we obtained 111 mRNAs (24 down-regulated and 87 up-regulated mRNAs). The predicted miRNAs and 251 lincRNAs were analysed using the 111 mRNAs, and 25 lincRNAs and 30 miRNAs were finally obtained. Subsequently, we used Cytoscape version 3.7.2 software to build a lincRNA-miRNA-mRNA ceRNA network (Fig. 5).

Fig. 5
figure5

The ceRNA network of lincRNA, miRNA and mRNA. Notes: Red hexagon denotes lincRNA, green triangle represents miRNA, and cyan ellipse represents mRNA

Functional enrichment analysis

A functional enrichment analysis was performed on the 111 mRNAs obtained previously. The biological processes of enrichment were mainly cell cycle, ossification, and sex differentiation (Fig. 6a). The cellular components were concentrated in the transcription factor, protein kinase, and serine/threonine protein kinase complexes (Fig. 6b). Enriched molecular function was mainly involved in DNA-binding transcription activator activity, RNA polymerase II-specific core promoter binding, and HMG box domain binding (Fig. 6c). KEGG pathway analysis showed that the genes were associated with cellular senescence, cell cycle, multiple cancers, miRNA in cancer, and TGF-beta signalling pathway (Fig. 6d).

Fig. 6
figure6

GO and KEGG analysis. a, The relationship between genes and GO terms of biological process. b, The relationship between genes and GO terms of cellular component. c, The relationship between genes and GO terms of molecular function. d, Chord plot indicates the relationship between genes and KEGG pathways. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes

Construction of prognostic models

Among the CGGA GBM samples, complete survival information was available for 236 samples, and all clinical information was available for 198 samples. Univariate Cox regression analysis was performed on the 111 mRNAs previously obtained, and 27 mRNAs related to survival time were obtained (p < 0.05) (Fig. 7a). Multivariate Cox analysis was performed using 27 mRNAs, and a Cox proportional hazards regression model of GBM patients containing 13 mRNAs was constructed (Fig. 7b). Based on the median risk score, all patients were divided into two groups (high-risk and low-risk groups). Figure 7c shows the survival status, survival time, and mRNA expression levels of patients. Survival analysis showed that patients in the low-risk group survived longer than those in the high-risk group (p < 0.001) (Fig. 7d). The univariate Cox regression analysis showed that recurrence (p = 0.004), isocitrate dehydrogenase (IDH) mutation (p = 0.009), 1p19q codeletion (p = 0.014), and risk score (p = 0.001) were predictors (Fig. 7e). Moreover, multivariate Cox regression analysis confirmed that recurrence (p = 0.002) and risk score (p < 0.001) were independent risk factors (Fig. 7f). The risk score and recurrence AUCs of the nomogram were 0.668 and 0.663, respectively (Fig. 7g).

Fig. 7
figure7

13-lincRNA prognosis model establishment and verification. a-b, Univariate Cox regression analysis and multivariate Cox regression analysis were used to construct prognostic models. c, Correlation between the prognostic signature and the overall survival of patients. The distribution of risk scores (upper), survival time (middle) and lincRNA expression levels (below). The black dotted lines represent the median risk score cut-off dividing patients into high- and low-risk groups. The red dots and lines represent the patients in high-risk groups. The green dots and lines represent the patients in low-risk groups. d, Kaplan-Meier survival curves of overall survival among risk stratification groups. e, Univariate Cox regression analyses of clinical factors associated with overall survival. f, Multivariate Cox regression analyses of clinical factors associated with overall survival. g, ROC curve related to clinical factors

CeRNA network in GBM

Using the 13 mRNAs obtained from the previous multivariate COX regression analysis conducted to screen the ceRNA network, a total of 14 miRNAs and 23 lincRNAs that interacted were found. We combined the previously selected 14 miRNAs and DEmiRNAs, screened out 6 overlapping miRNAs, and then screened the ceRNA network again. Finally, a ceRNA network containing 8 mRNAs, 6 miRNAs, and 18 lincRNAs was constructed (Fig. 8).

Fig. 8
figure8

Construction of a ceRNA network including 8 mRNAs, 6 miRNAs, and 18 lincRNAs. Notes: Red hexagon denotes lincRNA, green triangle represents miRNA, and cyan ellipse represents mRNA

GSEA reveals the close relationship between lincRNA and GBM

Using | log2FC | ≥2 as a condition to screen 18 kinds of lincRNA, we obtained seven of them: MALAT1, MEG3, NEAT1, MIR7-3HG, FAM95B1, EPB41L4A-AS1, and AC125494.1. Among them, extensive research has been conducted on MALAT1 and NEAT1. Thus, we selected MEG3, MIR7-3HG, FAM95B1, and EPB41L4A-AS1 in 248 CGGA samples for GSEA. As shown in Fig. 8, the GSEA results revealed that these lincRNAs were closely related to cancer. MEG3 showed potential effects on prostate cancer, colorectal cancer, and cancer pathways (Fig. 9a), MIR7-3HG demonstrated potential role in small cell lung cancer, prostate cancer, and cancer pathways (Fig. 9b), FAM95B1 was closely related to endometrial cancer, non-small cell lung cancer, renal cell carcinoma, and MTOR signalling pathway (Fig. 9c), while EPB41L4A-AS1 was associated with small cell lung cancer, prostate cancer, and JAK/STAT signalling pathway (Fig. 9d). AC125494.1 is not found in the CGGA data.

Fig. 9
figure9

Gene set enrichment analysis (GSEA). a, MEG3; b, MIR7-3HG; c, FAM95B1; d, EPB41L4A-AS1. P < 0.05

Discussion

GBM, which is the most common type of glioma, accounts for about 15% of all brain tumours. It has been reported that there are about three GBM patients per 100,000 people, and the average survival time of patients is only approximately 12–18 months. The 5-year and 10-year survival rates are approximately 5.5 and 2.9%, respectively [28, 29]. The low recovery rate and poor survival time of GBM patients may be related to the lack of efficient therapeutic targets. Therefore, it is of great significance to find new biomolecular markers and therapeutic targets of GBM. Currently, many scholars and research institutions focus on the role of non-coding RNAs in tumours, and several molecules have been significantly associated. Studies have shown that many lincRNAs can act as oncogenes or tumour suppressor genes in cancer. In recent years, ceRNA network-related research has gained intense attention. In a ceRNA network, lincRNAs can competitively bind miRNAs to regulate the expression of target mRNAs. Additionally, several studies have shown that ceRNA network is closely related to the occurrence and development of cancer and may be valuable in predicting the prognosis of patients.

In this study, we used WGCNA to find key modules related to GBM pathogenicity and combined differentially expressed genes for analysis. Based on the lincRNA-miRNA-mRNA interaction, a ceRNA network of GBM patients was preliminarily constructed, and 111 mRNAs in the network were functionally enriched. The results showed that these mRNAs have potential roles in cell cycle, various cancers, miRNA in cancer, DNA-binding transcriptional activator activity, TGF-β signalling pathway, transcription factor complex, RNA polymerase II specificity, etc. Further, univariate Cox regression and multivariate Cox regression analyses were used to construct a Cox proportional hazards regression model with 13 key genes: CDK6, FAM84B, FBLIM1, FJX1, GNB5, HOXA3, MLEC, PDXK, SOX11, SPRY4, TBPL1, TRIB2, and WEE1. We calculated the risk score of patients, and verified that the prognostic model is highly accurate. In addition, the nomogram based on the risk score has the highest AUC value. We combined DEmiRNA and 13 mRNAs to screen the ceRNA network. Finally, 8 mRNAs, 6 miRNAs, and 18 lincRNAs were used to construct a ceRNA network. We selected four lincRNAs, including MEG3, MIR7-3HG, FAM95B1, and EPB41L4A-AS1, and predicted their functions using the GSEA software. Our previous studies have shown that there are few studies on these four lincRNAs and GBM. These lincRNAs were shown to be closely related to a variety of cancers. Among them, MEG3 and MIR7-3HG are related to cancer pathways, while these two and EPB41L4A-AS1 all regulate the JAK/STAT signalling pathway.

MALAT1 and NEAT1 have been widely studied non-coding RNAs, and some studies have demonstrated that they have clear correlations with various cancers such as hepatocellular carcinoma and lung cancer [30,31,32,33]. It has also been reported that MALAT1 can be used as ceRNA of miR-199a to promote the expression of ZHX1, which in turn can regulate the proliferation of GBM cells [34]. The ceRNA effect between NEAT1 and miR-194-5p is related to the angiogenesis of glioma [35].

MEG3, also known as gene trap locus 2, is an imprinted gene located at 14q32 [36]. It is a new type of tumour suppressor that plays a role in various tumours, such as ovarian and bladder cancers [37, 38]. Studies have shown that MEG3 is under-expressed in gliomas, and its over-expression has a significant inhibitory effect on the proliferation and migration of glioma cells, while promoting its apoptosis and autophagy, and inhibiting the PI3K/AKT/mTOR signalling pathway [39, 40]. MEG3 can act as a ceRNA for miR-19a, thereby exerting an inhibitory effect on glioma [41]. The present study also predicted that MEG3 is potentially valuable in multiple tumours and cancer pathways.

According to previous reports, MIR7-3HG is associated with tumour progression and is highly expressed in endometrial cancer. Bioinformatics analysis by Wang et al. showed that its high expression in breast cancer is significantly related to the survival time of patients [42, 43]. It has also been reported that this gene can be used as ceRNA to up-regulate the expression of PEG10 by sponged miR-27a-3p, and thus plays a role in retinoblastoma [44]. EPB41L4A-AS1, located in the 5q22.2 region of the human genome, is an induced gene of p53 and pGC-1α, can regulate glycolysis and glutamine metabolism, and plays an important role in cancer metabolic reprogramming [45]. It is highly expressed in colorectal cancer tissues and is involved in the proliferation, invasion, and migration of colorectal cancer cells [46]. Additionally, FAM95B1 has been shown to be associated with thyroid cancer [47]. Limitations of our study include the lack of cell or animal experiments and its retrospective nature. Subsequently, we hope to verify our results by conducting cell-based experiments.

Conclusion

Conclusively, in this study, we constructed the lincRNA-miRNA-mRNA ceRNA network of GBM, which may be involved in its molecular regulation, and identified four lincRNAs with potential roles in the tumour. To our knowledge, three of these lincRNAs, namely MIR7-3HG, FAM95B1, and EPB41L4A-AS1, are novel potential therapeutic targets for GBM, as there are no previous related studies. Further, we created a prognostic model containing 13 genes, which may serve as a reliable prognostic indicator of GBM. Limitations of our study include the lack of cell or animal experiments and its retrospective nature. Subsequently, we hope to verify our results by conducting cell-based experiments. However, the present study is expected to facilitate in-depth understanding and study the molecular mechanism of GBM, and provide new insights into targeted therapy and prognosis of tumours.

Availability of data and materials

The data comes from TCGA, GTEx, GEO, CGGA databases, which are all public open platforms.

Abbreviations

LincRNA:

Long intergenic non-coding RNA

ceRNA:

Competitive endogenous RNA

GBM:

Glioblastoma

TCGA:

The cancer genome atlas

NCBI:

National center for biotechnology information

GEO:

Gene expression omnibus

CGGA:

Chinese glioma genome atlas

GTEx:

Genotype tissue expression

WGCNA:

weighted gene co-expression network analysis

GO:

Gene ontology

KEGG:

Kyoto encyclopedia of genes and genomes

GSEA:

Gene set enrichment analysis

References

  1. 1.

    Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, Ohgaki H, Wiestler OD, Kleihues P, Ellison DW. The 2016 world health organization classification of tumors of the central nervous system: a summary. Acta Neuropathol. 2016;131(6):803–20.

    PubMed  Article  Google Scholar 

  2. 2.

    Tanaka S, Louis DN, Curry WT, Batchelor TT, Dietrich J. Diagnostic and therapeutic avenues for glioblastoma: no longer a dead end? Nat Rev Clin Oncol. 2013;10(1):14–26.

    CAS  PubMed  Article  Google Scholar 

  3. 3.

    Djebali S, Davis CA, Merkel A, Dobin A, Lassmann T, Mortazavi A, Tanzer A, Lagarde J, Lin W, Schlesinger F, et al. Landscape of transcription in human cells. Nature. 2012;489(7414):101–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  4. 4.

    Batista Pedro J, Chang Howard Y. Long noncoding RNAs: cellular address codes in development and disease. Cell. 2013;152(6):1298–307.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  5. 5.

    Zhang HM, Zhu JK. Emerging roles of RNA processing factors in regulating long non-coding RNAs. RNA Biol. 2014;11(7):793–7.

    PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    Brunner AL, Beck AH, Edris B, Sweeney RT, Zhu SX, Li R, Montgomery K, Varma S, Gilks T, Guo X, et al. Transcriptional profiling of long non-coding RNAs and novel transcribed regions across a diverse panel of archived human cancers. Genome Biol. 2012;13(8):R75.

    PubMed  PubMed Central  Article  Google Scholar 

  7. 7.

    Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, Tilgner H, Guernec G, Martin D, Merkel A, Knowles DG, et al. The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. 2012;22(9):1775–89.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  8. 8.

    Ulitsky I. Bartel DP: lincRNAs: genomics, evolution, and mechanisms. Cell. 2013;154(1):26–46.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. 9.

    Tsai MC, Spitale RC, Chang HY. Long intergenic noncoding RNAs: new links in cancer progression. Cancer Res. 2011;71(1):3–7.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  10. 10.

    Pan Y, Li C, Chen J, Zhang K, Chu X, Wang R, Chen L. The emerging roles of Long noncoding RNA ROR (lincRNA-ROR) and its possible mechanisms in human cancers. Cell Physiol Biochem. 2016;40(1–2):219–29.

    CAS  PubMed  Article  Google Scholar 

  11. 11.

    Mello SS, Sinow C, Raj N, Mazur PK, Bieging-Rolett K, Broz DK, Imam JFC, Vogel H, Wood LD, Sage J, et al. Neat1 is a p53-inducible lincRNA essential for transformation suppression. Genes Dev. 2017;31(11):1095–108.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  12. 12.

    Zhou LN, Tian Y, Guo F, Yu B, Li J, Xu H, Su Z. LincRNA-p21 knockdown reversed tumor-associated macrophages function by promoting MDM2 to antagonize* p53 activation and alleviate breast cancer development. Cancer Immunol Immunother. 2020;69(5):835–46.

    CAS  PubMed  Article  Google Scholar 

  13. 13.

    Long J, Zhang L, Wan X, Lin J, Bai Y, Xu W, Xiong J, Zhao H. A four-gene-based prognostic model predicts overall survival in patients with hepatocellular carcinoma. J Cell Mol Med. 2018;22(12):5928–38.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Zhou JG, Zhao HT, Jin SH, Tian X, Ma H. Identification of a RNA-seq-based signature to improve prognostics for uterine sarcoma. Gynecol Oncol. 2019;155(3):499–507.

    CAS  PubMed  Article  Google Scholar 

  15. 15.

    Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi Pier P. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell. 2011;146(3):353–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  16. 16.

    Zhang M, Du X. Noncoding RNAs in gastric cancer: research progress and prospects. World J Gastroenterol. 2016;22(29):6610–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  17. 17.

    Wang L, Cho KB, Li Y, Tao G, Xie Z, Guo B. Long noncoding RNA (lncRNA)-mediated competing endogenous RNA networks provide novel potential biomarkers and therapeutic targets for colorectal cancer. Int J Mol Sci. 2019;20(22):5758.

  18. 18.

    Gao S, Wang P, Hua YQ, Xi H, Meng Z, Liu T, Chen Z, Liu L. ROR functions as a ceRNA to regulate Nanog expression by sponging miR-145 and predicts poor prognosis in pancreatic cancer. Oncotarget. 2015;7(2):1608–18.

    PubMed Central  Article  PubMed  Google Scholar 

  19. 19.

    Loewen G, Zhuo Y, Zhuang Y, Jayawickramarajah J, Shan B. lincRNA HOTAIR as a novel promoter of cancer progression. J Cancer Res Updates. 2014;3(3):134–40.

    Article  CAS  Google Scholar 

  20. 20.

    Wang H, Zhang H, Zeng J. Tan Y: ceRNA network analysis reveals prognostic markers for glioblastoma. Oncol Lett. 2019;17(6):5545–57.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Zhu X, Jiang L, Yang H, Chen T, Wu X, Lv K. Analyzing the lncRNA, miRNA, and mRNA-associated ceRNA networks to reveal potential prognostic biomarkers for glioblastoma multiforme. Cancer Cell Int. 2020;20:393.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. 22.

    Liu Z, Wang X, Yang G, Zhong C, Zhang R, Ye J, Zhong Y, Hu J, Ozal B, Zhao S. Construction of lncRNA-associated ceRNA networks to identify prognostic lncRNA biomarkers for glioblastoma. J Cell Biochem. 2020;121(7):3502–15.

    CAS  PubMed  Article  Google Scholar 

  23. 23.

    Carvalho BS, Irizarry RA. A framework for oligonucleotide microarray preprocessing. Bioinformatics. 2010;26(19):2363–7.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  24. 24.

    Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  25. 25.

    Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  26. 26.

    Yu GC, Wang LG, Han YY, He QY. ClusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  27. 27.

    Walter W, Sánchez-Cabo F, Ricote M. GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics. 2015;31(17):2912–4.

    CAS  PubMed  Article  Google Scholar 

  28. 28.

    Ostrom QT, Gittleman H, Stetson L, Virk SM, Barnholtz-Sloan JS. Epidemiology of Gliomas. In: Raizer J, Parsa A, editors. Current Understanding and Treatment of Gliomas, vol. 163. Cham: Springer International Publishing; 2015. p. 1–14.

    Google Scholar 

  29. 29.

    Ostrom QT, Gittleman H, Liao P, Vecchione-Koval T, Wolinsky Y, Kruchko C, Barnholtz-Sloan JS. CBTRUS Statistical Report: Primary brain and other central nervous system tumors diagnosed in the United States in 2010–2014. Neuro-Oncology. 2017;19(suppl_5):v1–v88.

    PubMed  PubMed Central  Article  Google Scholar 

  30. 30.

    Peng NF, He JR, Li JD, Huang H, Huang W, Liao Y, Zhu S. Long noncoding RNA MALAT1 inhibits the apoptosis and autophagy of hepatocellular carcinoma cell by targeting the microRNA-146a/PI3K/Akt/mTOR axis. Cancer Cell Int. 2020;20(1):165.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  31. 31.

    Zhu LY, Yang NH, Li CC, Liu G, Pan W, Li X. Long noncoding RNA NEAT1 promotes cell proliferation, migration, and invasion in hepatocellular carcinoma through interacting with miR-384. J Cell Biochem. 2019;120(2):1997–2006.

    CAS  Article  Google Scholar 

  32. 32.

    Song J, Su ZZ, Shen QM. Long non-coding RNA MALAT1 regulates proliferation, apoptosis, migration and invasion via miR-374b-5p/SRSF7 axis in non-small cell lung cancer. Eur Rev Med Pharmacol Sci. 2020;24(4):1853–62.

    CAS  PubMed  Google Scholar 

  33. 33.

    Zang FR, Rao YQ, Zhu XH, Wu Z, Jiang H. Shikonin suppresses NEAT1 and Akt signaling in treating paclitaxel-resistant non-small cell of lung cancer. Mol Med. 2020;26(1):28.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  34. 34.

    Liao KM, Lin YY, Gao WZ, Xiao Z, Medina R, Dmitriev P, Cui J, Zhuang Z, Zhao X, Qiu Y, et al. Blocking lncRNA MALAT1/miR-199a/ZHX1 axis inhibits glioblastoma proliferation and progression. Mol Ther Nucleic Acids. 2019;18:388–99.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  35. 35.

    Wang CL, Chen YX, Wang Y, Liu X, Liu Y, Li Y, Chen H, Fan C, Wu D, Yang J. Inhibition of COX-2, mPGES-1 and CYP4A by isoliquiritigenin blocks the angiogenic Akt signaling in glioma through ceRNA effect of miR-194-5p and lncRNA NEAT1. J Exp Clin Cancer Res. 2019;38(1):371.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  36. 36.

    Benetatos L, Vartholomatos G, Hatzimichael E. MEG3 imprinted gene contribution in tumorigenesis. Int J Cancer. 2011;129(4):773–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  37. 37.

    Buttarelli M, De Donato M, Raspaglio G, Babini G, Ciucci A, Martinelli E, Baccaro P, Pasciuto T, Fagotti A, Scambia G, et al. Clinical value of lncrna MEG3 in high-grade serous ovarian cancer. Cancers. 2020;12(4):966.

    CAS  PubMed Central  Article  Google Scholar 

  38. 38.

    Fan XR, Huang HF, Ji ZG, Mao QZ. Long non-coding RNA MEG3 functions as a competing endogenous RNA of miR-93 to regulate bladder cancer progression via PI3K/AKT/mTOR pathway. Trans Cancer Res. 2020;9(3):1678–88.

  39. 39.

    Zhao HK, Wang X, Feng XY, Li X, Pan L, Liu J, Wang F, Yuan Z, Yang L, Yu J, et al. Long non-coding RNA MEG3 regulates proliferation, apoptosis, and autophagy and is associated with prognosis in glioma. J Neuro-Oncol. 2018;140(2):281–8.

    CAS  Article  Google Scholar 

  40. 40.

    Xu DH, Chi GN, Zhao CH, Li DY. Long noncoding RNA MEG3 inhibits proliferation and migration but induces autophagy by regulation of Sirt7 and PI3K/AKT/mTOR pathway in glioma cells. J Cell Biochem. 2018.

  41. 41.

    Qin N, Tong GF, Sun LW, Xu XL. Long noncoding RNA MEG3 suppresses glioma cell proliferation, migration, and invasion by acting as a competing endogenous RNA of miR-19a. Oncol Res. 2017;25(9):1471–8.

    PubMed  Article  Google Scholar 

  42. 42.

    Dong PX, Xiong Y, Yue JM, Xu D, Ihira K, Konno Y, Kobayashi N, Todo Y, Watari H. Long noncoding RNA NEAT1 drives aggressive endometrial cancer progression via miR-361-regulated networks involving STAT3 and tumor microenvironment-related genes. J Exp Clin Cancer Res. 2019;38(1):295.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  43. 43.

    Wang JJ, Huang YQ, Song W, Li YF, Wang H, Wang WJ, Huang M. Comprehensive analysis of the lncRNA-associated competing endogenous RNA network in breast cancer. Oncol Rep. 2019;42(6):2572–82.

    CAS  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Ding FK, Jiang K, Sheng YJ, Li C, Zhu H. LncRNA MIR7-3HG executes a positive role in retinoblastoma progression via modulating miR-27a-3p/PEG10 axis. Exp Eye Res. 2020;193:107960.

    CAS  PubMed  Article  Google Scholar 

  45. 45.

    Liao MJ, Liao WJ, Xu NH, Li B, Liu F, Zhang S, Wang Y, Wang S, Zhu Y, Chen D, et al. LncRNA EPB41L4A-AS1 regulates glycolysis and glutaminolysis by mediating nucleolar translocation of HDAC2. EBioMedicine. 2019;41:200–13.

    PubMed  PubMed Central  Article  Google Scholar 

  46. 46.

    Bin J, Nie S, Tang Z, Kang A, Fu Z, Hu Y, Liao Q, Xiong W, Zhou Y, Tang Y, et al. Long noncoding RNA EPB41L4A-AS1 functions as an oncogene by regulating the rho/ROCK pathway in colorectal cancer. J Cell Physiol. 2020.

  47. 47.

    Li N, Cui MM, Yu P, Li Q. Correlations of lncRNAs with cervical lymph node metastasis and prognosis of papillary thyroid carcinoma. Onco Targets Ther. 2019;12:1269–78.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgments

Not applicable.

Funding

This study was funded by the Mudanjiang Medical University Graduate Innovation Fund of China (grant number 2019YJSCX-05MY), Hong Qi Scientific Research Fund Science and Technology Program (grant number 2019HQ-06). The funding body played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Affiliations

Authors

Contributions

GDL, DL and HG were involved in the concept and design of the study. GDL, JH, JL and CW participated in the data collection, analysis of data and preparation of the manuscript. GDL, GYL and SG carried out interpretation of data and drafted the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Haidong Gong.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Liu, G., Liu, D., Huang, J. et al. Comprehensive analysis of ceRNA network related to lincRNA in glioblastoma and prediction of clinical prognosis. BMC Cancer 21, 98 (2021). https://doi.org/10.1186/s12885-021-07817-5

Download citation

Keywords

  • Glioblastoma
  • ceRNA
  • lincRNA
  • Weighted gene co-expression network analysis
  • Prognostic model