Bioinformatic analysis of the potential molecular mechanism of PAK7 expression in glioblastoma.

The present study aimed to determine the potential molecular mechanisms underlying p21 (RAC1)‑activated kinase 7 (PAK7) expression in glioblastoma (GBM) and evaluate candidate prognosis biomarkers for GBM. Gene expression data from patients with GBM, including 144 tumor samples and 5 normal brain samples, were downloaded. Long non‑coding RNAs (lncRNAs) and messenger RNAs (mRNAs) were explored via re‑annotation. The differentially expressed genes (DEGs), including differentially expressed mRNAs and differentially expressed lncRNAs, were investigated and subjected to pathway analysis via gene set enrichment analysis. The miRNA‑lncRNA‑mRNA interaction [competing endogenous RNA (ceRNA)] network was investigated and survival analysis, including of overall survival (OS), was performed on lncRNAs/mRNAs to reveal prognostic markers for GBM. A total of 954 upregulated and 1,234 downregulated DEGs were investigated between GBM samples and control samples. These DEGs, including PAK7, were mainly enriched in pathways such as axon guidance. ceRNA network analysis revealed several outstanding ceRNA relationships, including miR‑185‑5p‑LINC00599‑PAK7. Moreover, paraneoplastic antigen Ma family member 5 (PNMA5) and somatostatin receptor 1 (SSTR1) were the two outstanding prognostic genes associated with OS. PAK7 may participate in the tumorigenesis of GBM by regulating axon guidance, and miR‑185‑5p may play an important role in GBM progression by sponging LINC00599 to prevent interactions with PAK7. PNMA5 and SSTR1 may serve as novel prognostic markers for GBM.

Introduction

Glioblastoma (GBM) represents 15% of all brain tumors (1). In the majority of cases, the length of survival following diagnosis is 12–15 months, with fewer than 3–5% of patients surviving for >5 years (2). In general, surgery after chemotherapy and radiation therapy are commonly used for GBM treatment (3). However, the cancer usually recurs, as there is no clear way to prevent the disease (4).

Improved understanding of the mechanisms underlying GBM at molecular and structural levels is valuable for clinical treatment (5). Bioinformatics may be effectively used to analyze GBM microarray data in order to provide a theoretical reference for further exploration of tumorigenic mechanisms and aid in searching for potential target genes (6). Based on bioinformatic studies, certain differentially expressed genes (DEGs), including p21 (RAC1)-activated kinase 7 (PAK7), were identified as potential therapeutic targets for glioma (7). As the longest member of the PAK family, PAK7 (renamed from PAK5) contains a total of eight protein-coding exons (8). PAK7, initially cloned as a brain-specific kinase, was hypothesized to be involved in neurite growth in neuronal cells (9,10). The downregulation of PAK7 expression was shown to inhibit glioma cell migration and invasion through competing endogenous RNA (ceRNA) networks such as the PAK7-early growth response 1-matrix metalloproteinase 2 signaling pathway (11). The levels of PAK7 are regulated by microRNAs (miRNAs/miRs) such as miR-129 and are associated with long non-coding RNAs (lncRNAs) in glioma cell progression (12). However, the potential molecular mechanisms underlying PAK7 activity in GBM remain unclear.

In the present study, gene expression data from GBM were downloaded and reanalyzed using bioinformatic tools. The DEGs, including differentially expressed lncRNAs (DElncRNAs) and differentially expressed mRNAs (DEMs) between GBM samples and normal samples were investigated and subjected to functional and pathway enrichment analyses. PAK7-associated miRNAs were investigated and a PAK7-associated ceRNA network was constructed. Furthermore, The Cancer Genome Atlas (TCGA) survival analysis was performed to confirm the effect of PAK7 expression on GBM. The study aimed to reveal the potential molecular mechanisms underling PAK7 expression in GBM and identify candidate prognostic biomarkers for GBM.

Materials and methods

Microarray data

The RNA-seq profile data and clinical data of TCGA-GBM were downloaded from the University of California Santa Cruz (UCSC) database (http://xena.ucsc.edu) based on the platform of Illumina HiSeq 2000 RNA Sequencing (Illumina, Inc.). A total of 144 tumor samples (GBM group) and five normal brain samples (control group) were enrolled in the present study.

Exploration of lncRNAs and mRNAs by re-annotation

An exon was defined as lncRNA or mRNA if its starting and ending positions (including the positive and negative chain) in the current RNAseq data were in accordance with the chromosome location annotation of lncRNA or mRNA (hg19, V25) in GENCODE database (https://www.gencodegenes.org/) (13), respectively.

Investigation of DEGs

Genes with low expression (expression value was 0 in >50% samples) were removed. Linear regression and empirical Bayesian methods (14,15) in Limma package of R (version: 3.40.2) (16) were used to explore DElncRNAs and DEMs. Multiple correction on the P-value of each DEG was performed using the Benjamini-Hochberg method (17) to obtain the adjusted P value. An adjusted P<0.05 and log|fold change| >1.5 were selected as the thresholds for DEG screening. A volcano plot was constructed based on the different expression values of DEMs and DElncRNAs.

Non-metric multidimensional scaling (NMDS) (18) was performed on DElncRNAs and DEMs using the vegan (version 2.5–2; http://cran.r-project.org/web/packages/vegan/) package (19) in R, followed by the construction of an NMDS plot. Analysis of similarities (20) was performed for sample similarity analysis. The distance between two points was defined as Bray-Curtis, and displacement tests were carried out. The threshold was set as P<0.05.

Gene set enrichment analysis (GSEA) for PAK7

Pearson's correlation coefficients for all mRNA-PAK7 pairs were calculated and pathway analysis (21) was performed on these mRNA-PAK7 pairs (according to the descending order of the correlation coefficient) via GSEA in the clusterProfiler package of R (version: 3.12.0) (22). A value of P<0.05 was considered as the threshold for current analysis.

miRNA prediction and co-expression investigation

Pearson's correlation coefficients for all DElncRNA-PAK7 and DEM-PAK7 relations were calculated and correlation tests were performed. DElncRNAs and DEMs with |r| >0.6 and P<0.05 were considered as key genes that were significantly associated with PAK7. The correlation between lncRNAs and mRNAs was calculated with the methods mentioned above. miRNA predictions for DEMs that positively correlated with mRNAs (including PAK7) was performed using miRWalk2.0 software (23,24) and combined with the prediction results from six databases, including miRWalk, miRanda (http://www.microrna.org/) (25), miRDB (http://www.mirdb.org/) (26), miRMap (http://www.mirmap.ezlab.org/) (27), RNA22 (https://www.rna-seqblog.com/rna22) (28) and TargetScan 7.1 (http://www.targetscan.org/vert_71/) (29). PAK7-DEMs-miRNAs that appeared in at least five databases were selected for further investigation.

Based on the miRNA-lncRNA relations in both starBase (30) (version 2.0; http://starbase.sysu.edu.cn/) and InCeDB (31) (http://gyanxet-beta.com/lncedb/) databases, the DElncRNAs that were positively related with PAK7 were predicted, and PAK7-DElncRNA-miRNA associations were investigated. The RNAs with same miRNA-binding site were defined as ceRNAs. Based on the PAK7-DEM-miRNA and PAK7-DElncRNA-miRNA interactions, the lncRNA-mRNA relationships regulated by the same miRNA were investigated. miR2Disease database (http://watson.compbio.iupui.edu:8080/miR2Disease/index.jsp) was used for miRNA investigation (search keyword: ‘Glioblastoma multiforme’) followed by lncRNA-mRNA association analysis. The miRNA-lncRNA-mRNA regulation network was constructed using lncRNA-mRNA and DEM-DElncRNA relationships obtained. Finally, lncRNAs and mRNAs regulated by the same miRNA (PAK7 was also regulated by this miRNA) were screened as key ceRNAs for PAK7.

Survival analysis

To reveal the prognostic value of the hub genes in patients with GBM, survival analysis, including overall survival (OS) and OS status, was performed. The samples in the data were divided into a high expression group (up group) and low expression group (down group) according to the mean and median values of all lncRNAs/mRNAs in the ceRNA network. The survival package (version 2.41–3; http://cran.r-project.org/web/views/Survival.html) in R was used for the log-rank test (32) in two groups. P<0.05 was considered statistically significant. Associations between lncRNAs/mRNAs and GBM prognosis were analyzed. The Kaplan-Meier (KM) (33) survival curve was constructed based on these selected lncRNAs/mRNAs.

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

A total of 8 tumor specimens were collected from patients with primary GBM (33–79 years old; 3 female and 5 male) who underwent neurosurgery at Fourth Affiliated Hospital of Harbin Medical University from February 2017 to July 2018, and 8 normal brain specimen samples (44–83 years; 4 female and 4 male) were obtained from patients who underwent cerebral contusion at Fourth Affiliated Hospital of Harbin Medical University from April 2017 to May 2018. Patients with a family history of cancer, and patients with chronic heart, lung, liver and kidney diseases were excluded. The study was approved by the Human Ethics Committee of Fourth Affiliated Hospital of Harbin Medical University and informed consent was obtained from all patients. RT-qPCR analysis was performed to detect the expression levels of several DElncRNAs (LINC00599) and DEMs [PAK7, PNMA family member 5 (PNMA5), and somatostatin receptor 1 (SSTR1)], as well as miR-185-5p. Briefly, total RNA was extracted from eight GBM samples and eight normal brain tissue samples using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.). cDNA was synthesized using PrimeScript™ RT Master Mix (Perfect Real Time) (Takara Biotechnology, Co., Ltd.). The amplification was performed SYBR Premix Ex Taq kit (cat. no. DRR041A; Takara Bio, Inc.) under the following conditions on an ABI ViiA 7 Real-Time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.): 50°C for 3 min, 95°C for 3 min, then 25 cycles of 95°C for 10 sec and 60°C for 30 sec, followed by dissociation curve analysis (60°C-95°C; 0.5°C increments for 10 sec). The primers used are shown in Table I. Relative expression levels were normalized to GAPDH and calculated with the 2−ΔΔCq method (34), and each sample was assayed in triplicate.

Table I.

Primers used for RT-qPCR.

Primer	Primer sequence (5′→3′)
PAK7-hF	GTAGTAGTTCCCCAGCGTGC
PAK7-hR	ACATGGTCTCAGAGGTCCAG
PNMA5-hF	TGCCCAGTCACATACCAGGAA
PNMA5-hR	CATACTTCGGCCCTCATCTTTC
SSTR1-hF	TGGCTATCCATTCCATTTGACC
SSTR1-hR	AGGACTGCATTGCTTGTCAGG
LINC00599-hF	CTAAACCCCTTGCCCCACAA
LINC00599-hR	AGGTTTTACAGGAGGGCAGC
hsa-miR-185-5p-R	GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTCAGGA
hsa-miR-185-5p-F	GGCTGGAGAGAAAGGCAGT
GAPDH-hF	TGACAACTTTGGTATCGTGGAAGG
GAPDH-hR	AGGCAGGGATGATGTTCTGGAGAG
U6-hF	CTCGCTTCGGCAGCACA
U6-hR	AACGCTTCACGAATTTGCGT
human-U6-RT	GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACAAAATATG
Downstream universal primer	GTGCAGGGTCCGAGGT

F, forward; R, reverse; RT, reverse transcription; qPCR, quantitative PCR; h, human; PAK7, p21 (RAC1) activated kinase 7; PNMA5, paraneoplastic antigen Ma family member 5; SSTR1, somatostatin receptor 1.

Data are presented as the mean ± SD, and statistical analysis was performed using GraphPad Prism 5.0 (GraphPad Software, Inc.). The difference between two groups was analyzed with a t-test, and a P<0.05 was considered to indicate a statistically significant difference.

Results

Investigation of DElncRNAs and DEMs

A total of 34 upregulated and 92 downregulated DElncRNAs were revealed between two groups. A total of 920 upregulated DEMs and 1,142 downregulated DEMs were observed. The volcanic diagrams for DElncRNAs and DEMs are shown in Fig. 1A and B. The expression values of DElncRNAs and DEMs were used to construct NMDS plots for these genes (Fig. 1C and D). The results showed that the green points (samples in GBM group) and red points (samples in control group) were distinctly separated, indicative of clear separation of samples based on the DEMs and DElncRNAs.

Figure 1.

Volcano plot and non-metric multidimensional scaling plot for the differentially expressed mRNAs and lncRNAs. (A) Volcano plot for differentially expressed lncRNAs. (B) Volcano plot for differentially expressed mRNAs. (C) Non-metric multidimensional scaling plot for differentially expressed lncRNAs. (D) Non-metric multidimensional scaling plot for differentially expressed mRNAs. For volcano plots, the red and blue points represent upregulated and downregulated genes, respectively. For non-metric multidimensional scaling plots, the red and green points represent the samples from the normal group and tumor group, respectively. lncRNA, long non-coding RNA.

Pathway enrichment analysis

GSEA pathway analysis was performed on the mRNAs in mRNA-PAK7 pairs. The results showed that the upregulated mRNAs were mainly enriched in 14 pathways, including neuroactive ligand receptor interactions, calcium signaling pathways and axon guidance (Fig. 2A-C). The downregulated mRNAs were mainly enriched in 35 pathways, including p53 signaling pathways, asthma and allograft rejection (Fig. 2D-F). According to the P-value, the detailed information for the top six upregulated and downregulated pathways is shown in Table II. The results of the pathway analysis indicated that PKA7 participated in axon guidance pathway (Fig. 3).

Figure 2.

Results of gene set enrichment analysis for outstanding pathways enriched by the differentially expressed mRNAs. (A) Neuroactive ligand receptor interaction pathway. (B) Calcium signaling pathway. (C) Axon guidance pathway. (D) p53 signaling pathway. (E) Asthma pathway. (F) Allograft rejection pathway.

Table II.

Top 6 pathways enriched by up- and downregulated mRNAs.

A, Upregulated mRNAs
Pathway ID	P-value (adjusted)	Enriched genes
KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION	3.32×10−2	GABRA3, CHRNA4, HRH3, CHRM1, GABRB3…
KEGG_CALCIUM_SIGNALING_PATHWAY	3.32×10−2	CACNA1E, CHRM1, HTR6, PLCB1, GNAL…
KEGG_AXON_GUIDANCE	3.32×10−2	PAK1, PAK3, PAK6, PAK7, L1CAM…
B, Downregulated mRNAs
Pathway ID	P-value	Enriched genes
KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION	3.32×10−2	GABRA3, CHRNA4, HRH3, CHRM1, GABRB3…
KEGG_CALCIUM_SIGNALING_PATHWAY	3.32×10−2	CACNA1E, CHRM1, HTR6, PLCB1, GNAL…
KEGG_AXON_GUIDANCE	3.32×10−2	PAK1, PAK3, PAK6, PAK7, L1CAM…

Figure 3.

Schematic of the axon guidance pathway enriched by PAK7. Copyright for use of this image was granted.

Co-expression network analysis

A total of 312 mRNAs and 22 lncRNAs showed positive correlation with PAK7. Following miRNA prediction, a total of 25,440 miRNA-lncRNA-mRNA regulatory interactions were reported from the screening of lncRNA-mRNA pairs regulated by the same miRNA. The selection of disease-related miRNAs revealed a total of 185 miRNA-lncRNA-mRNA relationships. Combined positive co-expression lncRNA-mRNA relationships revealed a total of 159 miRNA-lncRNA-mRNA regulatory associations. The ceRNA network was constructed based on network connectivity (Fig. 4). The results revealed 103 nodes (including 6 miRNAs, 92 mRNAs and 5 lncRNAs) and 281 interactions (including 115 miRNA-mRNA interactions, 7 miRNA-lncRNA interactions and 159 lncRNA-mRNA interactions), including miR-185-5p-LINC00599-PAK7 in the current ceRNA network.

Figure 4.

Competing endogenous RNA network constructed with miRNAs, lncRNAs, and mRNAs. Green diamonds represent lncRNAs, blue circles represents mRNAs (PAK7 is highlighted in yellow), red triangles represent miRNAs, red arrow represents miRNA-lncRNA regulatory relations, gray arrows represent miRNA-mRNA regulatory relations, green dotted lines represent mRNA-lncRNA co-expression relations. PAK7, p21 (RAC1)-activated kinase 7; miRNA/miR, microRNA; lncRNA, long non-coding RNA.

The mRNAs and lncRNAs regulated by the same miRNAs as PAK7 were investigated, followed by the construction of a core ceRNA network. As shown in Fig. 5, a total of 63 nodes (2 miRNAs, 58 mRNAs and 3 lncRNAs) and 171 interactions (62 miRNA-mRNA relations, 3 miRNA-lncRNA relations and 106 lncRNA-mRNA relations) were observed in the current core ceRNA network. These mRNAs and lncRNAs were regulated by hsa-mir-185-5p or hsa-mir-221-3p together with PAK7.

Figure 5.

Core competing endogenous RNA network associated with PAK7. Green diamonds represent lncRNAs, blue circles represent mRNAs (PAK7 is highlighted in yellow), red triangles represent miRNAs, red arrows represent miRNA-lncRNA regulatory relations, gray lines represent miRNA-mRNA regulatory relations, green dotted lines represent mRNA-lncRNA co-expression relations, red arrows represent relations with hsa-miR-185-5p, green arrows represent relations with hsa-miR-221-3p. PAK7, p21 (RAC1)-activated kinase 7; miRNA/miR, microRNA; lncRNA, long non-coding RNA.

Survival analysis for core ceRNAs

A total of 143 samples with OS prognosis information, and mRNA and lncRNA data were explored. PNMA5 was a unique mRNA associated with OS prognosis upon grouping of all mRNAs/lncRNAs based on the median values of expression (Fig. 6A). Hippocalcin-like 4 and SSTR1 were two outstanding mRNAs associated with OS prognosis following grouping of all mRNAs/lncRNAs based on the mean value of expression (Fig. 6B and C). KM analysis showed that PNMA5, Hippocalcin-like 4 and SSTR1 may be the factors associated with prognosis in patients with GBM.

Figure 6.

Kaplan-Meier survival curve analysis of the mRNAs potentially associated with glioblastoma multiforme prognosis. (A) Kaplan-Meier survival analysis for PNMA5. (B) Kaplan-Meier survival analysis for HPCAL4. (C) Kaplan-Meier survival analysis for SSTR1. PNMA5, paraneoplastic antigen Ma family member 5; HPCAL4, hippocalcin-like 4; SSTR1, somatostatin receptor 1.

RT-qPCR validation

The results of RT-qPCR showed that the mRNA levels of PAK7 and SSTR1 were significantly downregulated in GBM samples compared with control samples. Furthermore, LINC00599 expression was downregulated and miR-185-5p expression was upregulated in GBM samples (Fig. 7). These results were consistent with the bioinformatics analysis results. However, contrary to the bioinformatics analysis results, the mRNA expression of PNMA5 was significantly downregulated in GBM samples compared with control samples.

Figure 7.

Relative expression levels of GBM-associated RNAs. LINC00599, PAK7, PNMA5, SSTR1 and miR-185-5p were analyzed in GBM and control tissue samples. miR-185-5p, microRNA-185-5p; PAK7, p21 (RAC1)-activated kinase 7; PNMA5, paraneoplastic antigen Ma family member 5; SSTR1, somatostatin receptor 1. *P<0.05, ***P<0.001.

Discussion

GBM is the most common malignant brain tumor in adults (35). Targeted gene therapy such as PAK7 serves as a therapeutic strategy for glioma (11). However, the detailed mechanisms underlying the effects of PAK7 expression on GBM progression remain incompletely understood. In the present study, gene expression data of GBM were downloaded and reanalyzed using bioinformatic tools. The results revealed a total of 954 upregulated and 1,234 downregulated DEGs between GBM and control samples. These DEGs, including PAK7, were mainly enriched in pathways such as axon guidance. ceRNA network analysis revealed several outstanding ceRNA relationships, including miR-185-5p-LINC00599-PAK7. Moreover, PNMA5 and SSTR1 were two outstanding prognostic genes associated with OS in survival analysis.

Bioinformatic approaches have been widely used for the identification of genes that are associated with GBM. Wei et al (36) identified DEGs regulated by transcription factors in GBM by analyzing the GSE4290 dataset downloaded from the Gene Expression Omnibus (GEO) database. Zhang et al (37) selected GSE50161 from GEO database to identify potential oncogenes associated with GBM progression. Moreover, Liu et al (38) analyzed the glioma gene expression profile dataset GSE4290 for the identification of DEGs. These studies were performed based on datasets downloaded from GEO database. In particular, the expression profile data in GEO usually lack clinical data and may be unsuitable for survival analysis. In the present study, a dataset obtained from TCGA database that included clinical data was used, permitting survival analysis. In addition to DEM selection, the potential lncRNAs were also selected, and a regulatory network was constructed to explore the potential molecular mechanism of GBM.

Axon guidance (also called axon pathfinding) is a subfield of neural development related to the process by which neurons send out axons to reach the correct targets (39). Agrawal et al (40) reported the development of acute axonal polyneuropathy after the resection of GBM, indicative of the close relationship between axon guidance and GBM. Axon guidance is accomplished with the involvement of relatively few guidance molecules (41). Certain guidance molecules, such as PLEXIN-B2, may promote the tumorigenesis of GBM through axon guidance pathways (42). The dysregulation in gene expression implicated in the canonical axon guidance pathway influences the glioma cell cycle (43). Furthermore, PAK family members, including PAK1, PAK4 and PAK6, commonly participate in tumorigenesis as guide molecules (44–46). Hing et al (47) demonstrated that PAK was not only a critical regulator of axon guidance, but also served as a downstream effector of Dock in vivo. Axon guidance was shown to be involved in the PAK-small GTPase signaling in patients with schizophrenia (48). In the present study, PAK7 was one of the most outstanding DEGs enriched in axon guidance pathway. Thus, we speculate that PAK7 may participate in the tumorigenesis of GBM via the axon guidance pathway.

Studies have shown that miRNAs may suppress the process of tumorigenesis through the induction of cell differentiation and cycle arrest (49). It was shown that the negative expression of miR-185 in GBM may inhibit the proliferation of glioma cell lines (50). Thus, miR-185 is commonly used as a predictive biomarker for the prognosis of malignant glioma (51,52). miRNAs play important roles in tumor progression by targeting lncRNAs (53). Tran et al (54) reported that lncRNAs such as LINC00176 may regulate cell cycle and survival by titrating the tumor suppressor miR-185. Despite miRNA-lncRNA interactions, the lncRNA-mRNA regulatory relations are also vital for tumor progression (55). A comprehensive analysis on the functional miRNA-mRNA regulatory network revealed the association of miRNA signatures with glioma (56). As an mRNA, PAK7 is related with tumor progression by targeting of miRNAs, including miR-492 (57). Pan et al (58) showed that miR-106a-5p inhibited the migration and invasion of renal cell carcinoma cells by targeting PAK7. Neuronal activity regulates spine formation, at least in part, by increasing miRNA transcription, which in turn activates PAK actin remodeling pathway (59). However, the association between miR-185-PAK7 interactions and GBM progression is unclear. In the present study, ceRNA network analysis showed that the miR-185-5p-LINC00599-PAK7 axis was one of the regulatory relationships most strongly associated with GBM. Thus, it is hypothesized that the miR-185-5p may play an important role in GBM progression by sponging the PAK7-LINC00599 interaction.

The members of the PNMA family have been identified as onconeuronal antigens that are aberrantly expressed on cancer cells (60). Distinct functional domains of PNMA5 mediate protein-protein interactions, nuclear localization and apoptosis signaling in human cancer cells (60). As a novel tumor-associated gene, PNMA5 plays an important role in maintaining the malignant phenotype of hepatocellular carcinoma (61). However, studies on PNMA5 expression in glioma or GBM are very rare. SSTR1 acts at several sites to inhibit the release of hormones or secretory proteins (62). Kiviniemi et al (63) indicated that SSTR2 may serve as a prognostic marker in high-grade glioma. The activation of SSTR1 was shown to induce glioma growth arrest in vitro and in vivo (64). However, whether SSTR1 can be used as a prognostic marker in GBM is unclear. In the present study, the survival analysis results showed that PNMA5 and SSTR1 were the two most notable DEGs associated with overall survival in disease. This result indicated that PNMA5 and SSTR1 may be used as novel prognostic markers for GBM. However, certain limitations of this study include the small sample size and lack of verification. Thus, a larger sample size with more extensive verification analysis is warranted in future investigations.

In conclusion, the present study suggested that PAK7 expression was downregulated in GBM, and that this gene may participate in the tumorigenesis of GBM through the ceRNA pathway of miR-185-5p-LINC00599-PAK7. Although the expression levels of miR-185-5p, LINC00599 and PAK7 in GBM tissue were detected, the regulatory relation of this ceRNA needs to be validated in future studies.