RNA sequencing and clinical data of GBM patients were obtained from the CGGA database (http://www.cgga.org.cn). The gene expression levels were measured in terms of fragments per kilobase of transcript per million mapped reads. Various clinical data were also downloaded from the database, which included the patient gender, age, The Cancer Genome Atlas (TCGA) subtype, overall survival (OS), radiotherapy and chemotherapy details, and the mutation status of the genes isocitrate dehydrogenase (IDH), tumour protein p53 (TP53), EGFR, ATRX and enhancer of zeste homolog 2 (EZH2) (11,12). All cases with pathological diagnosis of GBM were included in the analysis. The exclusion criteria applied in the present study were as follows: i) Histologic confirmation of the diagnosis of any brain tumour type other than primary GBM; ii) history of radiotherapy or chemotherapy prior to histologic diagnosis; iii) patients with missing follow-up records; and iv) missing mutation information for the five aforementioned genes. According to these criteria, a total of 88 GBM samples were selected for inclusion in the current study (Table I). Furthermore, another gene expression and clinical dataset of GBM patients was downloaded from the TCGA database (https://cancergenome.nih.gov/), from which 162 GBM cases were selected. The dataset obtained from TCGA was analysed, and served as the validation dataset.

First, the expression data profiles from the CGGA were tested to confirm that they were suitable for the analysis. The standard deviation value for all samples of each gene was calculated, and the top 5,000 genes with the lowest standard deviation values were selected for subsequent analysis. Next, the co-expression network of genes was constructed using the WGCNA package in R (8). To calculate the scale-free topology fitting index r² that corresponded to different soft-thresholding parameter β values, functional pickSoftThreshold was used, and the β value was selected if r² reached 0.9. The soft-thresholding power β value was then set to 6, and the minModuleSize was set to 30. Subsequently, the gene expression profile was transformed into an adjacency matrix and a topological overlap matrix (TOM), which was defined as the sum of adjacency between the gene and all other genes for network generation. Next, the corresponding dissimilarity of TOM (dissTOM) was calculated, and dissTOM-based hierarchical clustering was used to produce a hierarchical clustering dendrogram of genes. Modules of clustered genes were then generated using the Dynamic Tree Cut algorithm. The module eigengene was calculated using the function moduleEigengenes, and a number of modules were merged according to a cut-off line for the module dendrogram. The interactions (correlations) of each module were analysed and visualised by heatmaps. To identify modules associated with patient characteristics, the Pearson's correlation test was used to evaluate the correlation of module eigengenes with the clinical traits, OS and mutation status, and correlations with P-values of <0.05 were considered to be statistically significant.

The functional enrichment of the genes of the identified module was assessed based on GO terms (13) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway (14) annotations. GO term analyses were performed using the DAVID database (https://david.ncifcrf.gov/) and Panther database (http://www.pantherdb.org) (15), which are essential tools for the success of high-throughput gene function analysis. Pathway analysis was also conducted using multiple online databases, including the DAVID database, KEGG pathway database (http://www.genome.jp/kegg) and STRING online database (http://string-db.org) (16). P-values of <0.05 were considered to denote statistically significant differences in GO term enrichment and KEGG pathway analyses, and the false discovery rate was utilised to correct the P-values.

To identify the gene-encoded proteins and construct the PPI network of the identified module, the genes were mapped to the STRING database. The results obtained from this database were then imported into Cytoscape software (version 3.6.0; https://cyto-scape.org/) to analyse the interactional associations among the gene-encoding proteins and their degrees in GBM (17). In addition, significant genes from the PPI network complex were selected according to their degree of importance. The corresponding proteins may be the core proteins or key candidate genes that have significant physiological regulatory functions.

To confirm the reliability of the identified genes from the CGGA data, GBM data from TCGA were then used to perform validation with the GEPIA database (http://gepia.cancer-pku.cn) (18). Through this database, the expression levels of all genes of interest in GBM and other tumours can be obtained. Furthermore, Kaplan-Meier curves were generated based on the GEPIA database. The OS was estimated using the log-rank test, and P<0.05 was considered to denote statistically significant data.

In TCGA data, mRNAs having a Spearman's correlation with lncRNA of >0.4 were considered to be lncRNA-correlated mRNAs. These were then analysed by KEGG pathway enrichment analysis. A P-value of <0.05 was applied to identify the significant pathways.

The genes nearby lncRNAs were analysed by genomic region enrichment of annotations tool (GREAT version 3.0.0) (19). The potential targets of lncRNA were predicted by searching the miRDB database (http://www.mirdb.org).

To detect and explore the possible biological function of the key survival-associated genes, WGCNA was performed based on the mRNA and lncRNA profiles derived from the CGGA database. According to the exclusion criteria mentioned earlier, RNA sequencing results and the clinical data of 88 GBM samples were downloaded from the CGGA database. For module detection, 5,000 coding and non-coding RNAs were selected for further analysis from the original 21,000 genes according to the standard deviation values. One outlier sample was removed from the sample network. The TCGA subtype, gender, age, OS, radiotherapy and chemotherapy information, and the mutation status of IDH, TP53, EGFR, ATRX and EZH2 were defined as clinical traits (Fig. 1A). Analysis of the network topology was first performed for various soft-thresholding power β values to determine the relative balanced scale independence and mean connectivity of the weighted gene co-expression network. As shown in Fig. 1B, power 6 was the lowest power at which the scale-free topology fitting indices r² reached 0.90; thus, this power was selected in order to produce a hierarchical clustering tree (dendrogram) of the 5,000 genes. In total, 21 modules were identified by hierarchical clustering and dynamic branch cutting, and each module was assigned a unique colour as an identifier (Fig. 1C and D). The largest module contained 826 genes, while the smallest contained 58 genes. The grey module represented a gene set that was not assigned to any of the modules.

To explore the survival significance of the selected modules, correlations between the OS and module eigengenes were analysed. It was observed that three modules (black, brown and light cyan) were positively correlated with OS (P<0.05). Among them, the black module had the lowest P-value. Another five modules (turquoise, blue, purple, tan and grey) were negatively correlated with OS (Fig. 2; P<0.05), with the purple module exhibiting the lowest P-value. According to the previous analysis of the survival significance of modules, the black and purple modules were selected for further analysis with OS, since they had the lowest P-values. Limited progress in targeted therapy for GBM has been made in recent years, and the majority of the therapeutic targets for this disease are oncogenes (20); thus, the purple module was further analysed in the present study, since the genes in purple module are considered to be oncogenes. In total, 195 genes, including 193 protein-coding and 2 non-coding genes, were identified in the purple module. The locations of certain of these identified genes on human chromosomes are displayed in the Circos plot in Fig. 3.

The functions and pathway enrichment of the candidate genes in GBM identified from the purple module were analysed using multiple online databases, including the DAVID, Panther, KEGG pathway and STRING databases. The GO functional enrichments of genes were determined with the DAVID and Panther databases, with a P-value of <0.05 indicating statistical significance of the data. The protein-coding genes in the purple module were mapped to the GO database to determine their potential functions. GO terms were divided into three functional groups, including biological processes (BP), cell composition (CC) and molecular function (MF). The enriched GO terms for the genes are presented in Fig. 4. For the candidate genes in the purple module, the top three enriched GO terms in each category were as follows: Extracellular structure organisation, extracellular matrix (ECM) organisation and multicellular organismal macromolecule metabolic process in the BP category; ECM, extracellular region and extracellular space in the CC category; and collagen binding, ECM structural constituent and growth factor binding in the MF category (Fig. 4A and B). These results revealed that the majority of the genes were significantly enriched in extracellular structure, binding, cell parts and cell growth.

Subsequently, functional and signalling pathway enrichment of genes in the purple module was performed using the online databases DAVID, STRING and KEGG. The top enriched KEGG pathways for the candidate genes included ECM organisation, collagen formation, integrin cell surface interactions, degradation of the ECM, collagen biosynthesis and modifying enzymes, ECM proteoglycans, assembly of collagen fibrils and other multimeric structures, ECM-receptor interaction, collagen degradation and signalling by receptor tyrosine kinases (Fig. 5). The majority of the identified pathways were cancer-associated signalling pathways.

Data from the STRING database revealed that several of the genes interacted with each other. In total, 113 of the 193 candidate protein-coding genes were filtered to form the PPI network complex (Fig. 6). The network contained 113 nodes and 648 edges. Among the 113 nodes, the most significant 19 hub node genes were identified using a degree of ≥10 as the filtering criterion. These genes were COL1A1, COL1A2, DCN, COL3A1, FN1, MMP9, COL6A1, COL5A2, COL6A2, ITGA5, FBN1, CTGF, ITGA4, PDGFRB, LUM, PPIB, MMP14, ITGA11 and COPB1 (Table II). Among these genes, COL1A1 exhibited the highest node degree (degree =30), and FN1 shares the highest degree with 33 validated coding genes.

To confirm the reliability of the 193 survival-associated coding genes identified from the CGGA, GBM datasets were also downloaded from TCGA database (including 162 GBM samples), and the RNA sequencing data and survival information of these datasets were subjected to Kaplan-Meier survival analysis. The results revealed significantly different OS between the high and low expression groups for 33 genes in the TCGA GBM, including ADAM12, B4GALT7, CD248, CHPF2, COL6A1, COL6A2, CYGB, DCBLD2, DERL2, DUSP6, EFNB2, EMILIN1, EPHA2, FAP, FBLN1, FN1, FZD1, HOXB2, HSP90B1, IGFBP4, LAMB1, LOXL1, MMP11, NID2, P4HB, PCOLCE, PDIA3, PLOD1, PMAIP1, PPIB, RARRES1, TBL2 and THY1 (P<0.05; Fig. 7). Therefore, the expression levels of these 33 genes may be used as predictors of OS in GBM patients, indicating that they may be candidate genes involved in GBM that deserve further investigation.

According to the co-expression network, two lncRNAs, namely LOC541471 and LOC284494, were identified as survival-associated key non-coding genes in GBM. Next, the genes nearby these two lncRNAs were analysed. The results demonstrated that 8 genes were located near the transcription start site (Fig. 8A and B). To confirm the reliability of the identified lncRNAs, TCGA GBM data were subsequently used to perform validation by GEPIA. As shown in Fig. 8C, LOC541471 was significantly overexpressed in the GBM data-sets obtained from TCGA. Next, the expression levels of these two lncRNAs in other types of cancer were explored using TCGA. It was observed that LOC541471 was overexpressed in 17 other cancer types (Fig. 8D), while LOC284494 was overexpressed in 3 other cancer types, as compared with their corresponding normal tissues (data not shown). Subsequently, survival analysis was performed in the TCGA GBM data, and a significant association with survival was observed only for LOC541471 (Fig. 8E). Therefore, LOC541471 was regarded as a core lncRNA in the network and warrants further investigation.

As mentioned earlier, LOC541471 was found to be the core lncRNA in the co-expression network. To understand how this lncRNA is involved in GBM, a deeper insight into the GBM expression data from TCGA was required. Spearman's correlation of LOC541471 was calculated with ~24,300 coding genes in 156 GBM patients. It was observed that pleckstrin-2 (PLEK2) is the mRNA exhibiting the highest correlation index (r=0.68) with lncRNA LOC541471, and this gene is located on human chromosome 3p24.1 (Fig. 9A and B). According to the TCGA data, PLEK2 is significantly overexpressed in GBM tissues as compared with the normal tissue (Fig. 9C). Next, a total of 436 mRNAs that had a Spearman's correlation with LOC541471 of >0.4 were further examined by KEGG pathway analysis. Pathways were considered as enriched at a P-value of <0.05 (Fig. 9D), and the most enriched pathway was oxidative phosphorylation (OXPHOS). Finally, the expression data of GBM samples obtained from TCGA and CGGA were analysed, and the expression of LOC541471 in IDH1 wild-type was observed to be much higher compared with that in the IDH1 mutant group (Fig. 10).

Cytoplasmic lncRNAs can act as competing endogenous RNAs to modulate the functions of miRNAs (21). By searching the miRDB database (22), it was observed that LOC541471 was located in the cytoplasm, and 37 possible targets miRNAs were identified (Table III). The top 5 target miRNAs were hsa-miR-548t-3p, hsa-miR-548ap-3p, hsa-miR-548aa, hsa-miR-4288 and hsa-miR-3138. Taken together, these results predicted the potential mechanism of LOC541471; however, further studies are required to demonstrate this additional mechanism of LOC541471 in GBM.