Identification of Prognostic Biomarkers for Glioblastoma Based on Transcriptome and Proteome Association Analysis.

OBJECTIVE: Glioblastoma multiforme (GBM) is the most malignant primary brain tumor in adults. This study aimed to identify significant prognostic biomarkers related to GBM.
METHODS: We collected 3 GBM and 3 healthy human brain samples for transcriptome and proteomic sequencing analysis. Differentially expressed genes (DEGs) between GBM and control samples were identified using the edge R package in R. Functional enrichment analyses, prediction of long noncoding RNA target genes, and protein-protein interaction network analyses were performed. Subsequently, transcriptomic and proteomic association analyses, validation using The Cancer Genome Atlas (TCGA) database, and survival and prognostic analyses were conducted. Then the hub genes directly related to GBM were screened. Finally, the expression of key genes was verified by quantitative polymerase chain reaction (qPCR).
RESULTS: Totally, 1140 transcripts and 503 proteins were significantly up- or down-regulated. A total of 25 genes were upregulated and 62 were downregulated at both the transcriptome and proteome levels. Results from TCGA database showed that 84 of these 87 genes matched with transcriptome sequencing results. A Cox regression analysis suggested that Fibronectin 1(FN1) was a prognostic risk factor. The qPCR results showed that FN1 was significantly upregulated in GBM samples.
CONCLUSIONS: FN1 may play a role in GBM progression through ECM-receptor interaction and PI3K-Akt signaling pathways. FN1 may be considered as a prognostic biomarkers related to GBM.

Introduction

Glioma is one of the most common tumors originating in the central nervous system. Although the treatment of glioma has made great progress in the past few decades, the clinical outcome of most patients is still poor, especially for glioblastoma (GBM), with a median survival of only 14.6 months, after timely and effective treatment, only 5% to 10% of patients can survive for 2 years. In recent years, surgery combined with radiotherapy and tinidazolamide chemotherapy have become the standard treatment for glioblastoma, but even if they receive the same treatment, the survival of patients is significantly different. Survival is affected by many factors, including treatment strategy, patient’s physical condition, and tumor characteristics. Biomarkers related to these factors are often the key to evaluating the prognosis of patients with glioma. Targeted therapy for related markers will become a major breakthrough in the treatment of glioma in the future.

Accordingly, several previous studies have focused on the identification of potential biomarkers of GBM. Epidermal growth factor receptor and transforming and vascular endothelial growth factors have been identified as prognostic and diagnostic markers of GBM. Wei et al screened out 5 glioblastoma transcription factors and their downstream DEGs through bioinformatics analysis, and speculated that they may be potential biomarkers for clinical treatment of glioblastoma. Long et al predicted and analyzed multiple key genes and pathways of glioblastoma based on bioinformatics. Although some key genes related to the progression of GBM have been identified, most are based on the analysis of microarray data downloaded from public databases. There is currently no clear understanding of the exact molecular mechanism underlying GBM progression.

This study used high-throughput transcriptomics and proteomics methods to determine the prognostic biomarkers associated with GBM. We performed transcriptome and proteome sequencing on GBM and healthy control human brain samples, identified DEGs at the transcriptome level and differentially expressed proteins (DEPs) at the proteome level. The differentially expressed genes/proteins were analyzed by correlation and the consistency of the transcriptome and proteome levels were observed. The expression of related genes was further verified through the TCGA database, and the survival prognosis analysis and experimental verification of the genes with the same differential expression were carried out. Through a comprehensive analysis of these transcriptome and proteome data, it is expected to reveal the pathogenesis of GBM and provide a basis for improving the prognosis of GBM.

Material and Methods

Samples and Sequencing

Three GBM (T1, T2, and T5) and 3 healthy control (T3, T4, and T6) human brain samples were collected from our hospital. Patient demographic information is summarized in Table 1. Gene sequencing of the samples was performed and the corresponding data were uploaded to the National Center for Biotechnology Information database (SRP108388). Quality control was completed using the FASTX-Toolkit version 0.0.13 and PRINSEQ-Lite version 0.20.4. The length of the lower base call accuracy (Q20) and raw and clean reads were calculated to visually demonstrate the validity and reliability of sequencing results. After preprocessing, clean reads were mapped to human reference sequence hg19 using TopHat version 2.0.8, allowing for a read mismatch of 2, a read gap length of 2, and a min-anchor-length of 8.

Table 1.

The Demographic Information on the Patients.

Samples	Male/female	Age	Body weight/kg	Pathogenetic locations	Tumor size/cm3
GBM group
T1	Female	45	50	Left frontal	5*6*8
T2	Male	38	54	Right frontal	5*6*9
T5	Male	27	76	Left temporal	5*6*7
Control group
T3	Male	51	65	Right temporopartial	/
T4	Male	70	75	Right temporal	/
T6	Male	27	85	Left temporal	/

Abbreviation: GBM, glioblastoma.

All the procedures in our study were approved by the institutional ethics committee of The First Affiliated Hospital of Harbin Medical University and the methods were carried out in accordance with the approved guidelines. All the patients have been informed and signed informed consent before the experiments.

Analyses of DEGs

Based on the annotation of known genes and lncRNAs provided by GENCODE version 24 (mapped to GRCh37), FPKM value and read count of genes were obtained using StringTie version 1.2.2. The edgeR package in R was used to obtain DEGs with |log2FC| >1 and an adjusted P value of < 0.05.

Functional Enrichment Analysis

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for DEGs were performed using the clusterProfiler package in R.

Prediction of lncRNA Target Genes

Weighted Gene Coexpression Network Analysis (WGCNA) was performed using the WGCNA package. The association between every module and disease was obtained using ANNOVAR with a cutoff P value of <0.05. We assumed that the microRNA was the target gene of the lncRNA in the same network module.

PPI Network Analysis

PPI network analysis for DEGs was conducted using STRING, and the PPI network was visualized using Cytoscape. The MCODE plugin for Cytoscape was used to obtain a subnetwork module with a degree cutoff of 5, node score cutoff of 0.2, and a K-core value of 2.

Mass Spectrometry (MS) Data Processing

The AB SCIEX 5600 Plus mass spectrometer was used to analyze the 6 brain samples. The identification of peptides and proteins from MS/MS data was conducted using ProteinPilot version 4.5 software. An unused score of ≥1.3 (confidence level >95%) and proteins that contained at least one unique peptide were regarded as credible proteins; quantitative results with a confidence level <95% were filtered out. Control samples were marked with numbers 115, 116, and 121, whereas GBM samples were marked with numbers 117-119.

For the experimental design with biological and technical repetition, the mean values of pairwise comparison results were normalized and regarded as the fold-change of samples. Subsequently, the minimum P value of a Student’s t-test for pairwise comparison results was regarded as a significant difference between samples. Lastly, significantly differentially expressed proteins were screened according to fold-change and P values. Proteins with a fold-change of ≥1.5 (upregulated, ≥1.5 and downregulated, ≤0.67) and a P value of ≤0.05 were considered significantly differentially expressed.

Transcriptomic and Proteomic Association Analyses

Association analyses for DEGs and differentially expressed proteins were performed, and the consistency between transcriptome and proteome levels was observed. The GO biological process for coexpressed genes and proteins was analyzed using the Clue GO plugin for Cytoscape with a P value of <0.05 and a maximum level of <3.

Validation Using TCGA Database

The GBM RNASeq version 2 data downloaded from TCGA database was used to screen the expression of genes associated with the transcriptome and proteome. The read count of every gene was downloaded and DEGs for GBMs between controls were identified using the edge R package. |log2FC| >1 and false discovery rate <0.05 were regarded as cutoff criteria and the difference in the expression of DEGs in the transcriptome and proteome was observed.

Survival and Prognostic Analyses

The effects of DEGs on patient prognosis were assessed according to GBM survival data provided by TCGA database. First, one-way analysis of variance (ANOVA) was conducted for each gene to obtain potential prognostic genes. Subsequently, multivariate ANOVA was conducted for all potential prognostic genes; Akaike information criteria were introduced to optimize included variables. Lastly, multivariate ANOVA was conducted for optimized variables to further select for prognostic genes. The optimized prognostic genes were classified into 2 groups based on their expression level in patients and a Kaplan-Meier curve was constructed to compare survival between the high and low expression groups.

Experimental Verification

Three GBM samples and 3 control (paracancerous) samples were used for detection. Based on our previous studies, we measured the expression of Fibronectin 1(FN1). The primers for FN1 were as follows: FN1-hF, 5′-CGGTGGCTGTCAGTCAAAG-3′; FN1-hR, 5′-AAACCTCGGCTTCCTCCATAA-3′; human-GAPDHF, 5′-AGAAGGCTGGGGCTCATTTG-3′; human-GAPDHR, 5′-AGGGGCCATCCACAGTCTTC-3′.

The total RNA was extracted from samples using RNAiso Plus (Trizol) (TaKaRa#9109) according to the instructions. The reverse transcription reaction mixture was prepared as follows: 4-µL 5x primeScript RT Master MIX (perfect Real Time, TAKARA#RR036A) and 0.5µg Total RNA were mixed, and RNase free water was added up to 20 μL. The reaction conditions were 37°C for 60 min and then 85°C for 5 s. Then, the qPCR reaction solution was prepared with the following components: 10 μL SYBR Premix EX Taq (2x, Thermo#4367659), 1 μL forward primer 10μM, 1 μL reverse primer 10 μM, and 8 μL cDNA. The qPCR reaction was performed using the following steps: 50°C for 3 min, 40 cycles of 95°C for 3 min, 95 for 10 s, and 60°C for 30 s. Finally, melt curve analysis was carried out in 60°C-95°C, using increments of 0.5°C per 10 s.

All results are presented as the mean ± SEM and presented in tables. SPSS version 22.0 software (IBM Corp., Armonk, NY, USA) was used for the statistical analyses, and GraphPad Prism 5 (GraphPad Software, San Diego, CA, USA) was used for mapping. P < 0.05 and P < 0.01 were considered significant and extremely significant differences, respectively.

Results

Compared with control samples, 1140 DEGs (332 upregulated and 808 downregulated) were identified in GBM samples, 102 of which were lncRNAs (17 upregulated and 85 downregulated). The heatmaps for DEGs and differentially expressed lncRNAs are shown in Figure 1.

Figure 1.

Heatmaps for differentially expressed genes (DEGs) (A) and differentially expressed long noncoding RNAs (lncRNAs) (B).

As shown in Figure 2, upregulated DEGs (e.g., FN1) were mainly enriched in ECM–receptor interaction pathways, the cell cycle, and PI3K–Akt signaling. Downregulated DEGs were mainly enriched in glutamatergic synapses, calcium signaling, and neuroactive ligand–receptor interactions.

Figure 2.

Functional enrichment analysis of upregulated (A) and downregulated (B) differentially expressed genes (DEGs).

In total, 4323 edges and 773 nodes (231 upregulated and 542 downregulated) were included in the PPI network. The top 10 degrees of nodes in the PPI network included tumor protein 53 (degree = 56), cyclin dependent kinase 1 (degree = 53), vascular endothelial growth factor A (degree = 49), cyclin B1 (degree = 46), kinesin family member 11 (degree = 44), baculoviral IAP repeat containing 5 (degree = 44), cell division cycle 20 (degree = 42), polo-like kinase 1 (degree = 41), MYC proto-oncogene (degree = 41), cyclin dependent kinase 2 (degree = 41), and aurora kinase A (degree = 41).

Network Module Analysis

In total, 14 network modules were obtained from the PPI network. The top 4 modules are displayed graphically in Figure 3A. KEGG pathway enrichment analyzes showed that Module 1 was significantly enriched in the cell cycle, oocyte meiosis, and neuroactive ligand–receptor interactions; Modules 2 and 3 were significantly enriched in neuroactive ligand–receptor interactions; and Module 4 was significantly enriched in protein digestion and absorption, ECM–receptor interaction pathways, and focal adhesion (Figure 3B).

Figure 3.

Network module analysis results of the top 4 modules (A), and their Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyzes (B).

WGCNA

The DEGs were divided into WGCNA modules, and each module was set to contain at least 10 genes. All differentially expressed genes were divided into 11 modules, and the correlation between the modules and tumor was obtained by single factor ANOVA analysis. The results regarding the association between every module and disease suggested that the lightcyan and black modules were significantly associated with GBM (Table 2). There were 560 genes, including 49 lncRNAs, in the lightcyan module, and 162 genes, including 13 lncRNAs, in the black module. KEGG pathway enrichment analyzes for genes excluding lncRNAs indicated that the black module was significantly enriched in long-term potentiation, whereas the lightcyan module was significantly enriched in neuroactive ligand–receptor interactions and glutamatergic synapses (Table 3). The network for the top 50 genes in modules is presented in Figure 4.

Table 2.

Association Results of WGCNA Network Module and Disease State.a

Module	MEmagenta	MElightcyan	MEblack	MElightgreen	MEdarkred	MEblue	MEdarkgrey	MEdarkturquoise	MEgrey
Cor	−0.7	−0.93	−0.88	−0.56	0.63	0.66	0.67	0.8	−0.3
P vlaue	0.124462901	0.007922162	0.020543815	0.243587746	0.184230176	0.1573967	0.144167288	0.055461844	0.559032811

a The boldface values suggest that the lightcyan and black modules were significantly associated with GBM.

Table 3.

Enrichment Results of KEGG Metabolic Pathways of the 2 Modules.

Module	ID	Description	P value	Count
Black	hsa04720	Long-term potentiation	8.11E-06	6
	hsa04724	Glutamatergic synapse	1.68E-05	7
	hsa04072	Phospholipase D signaling pathway	6.34E-04	6
Lightcyan	hsa04080	Neuroactive ligand-receptor interaction	3.47E-10	28
	hsa04724	Glutamatergic synapse	4.77E-10	18
	hsa04020	Calcium signaling pathway	1.04E-09	22
	hsa04727	GABAergic synapse	4.01E-08	14
	hsa05033	Nicotine addiction	4.22E-08	10
	hsa05031	Amphetamine addiction	1.01E-06	11
	hsa05032	Morphine addiction	3.07E-06	12
	hsa05030	Cocaine addiction	3.36E-06	9
	hsa04721	Synaptic vesicle cycle	3.81E-06	10
	hsa04728	Dopaminergic synapse	5.31E-06	14
	hsa04713	Circadian entrainment	5.43E-06	12
	hsa04024	cAMP signaling pathway	1.15E-05	17
	hsa04725	Cholinergic synapse	2.68E-05	12
	hsa04921	Oxytocin signaling pathway	3.47E-05	14
	hsa04360	Axon guidance	3.92E-05	15
	hsa04911	Insulin secretion	5.72E-05	10
	hsa04014	Ras signaling pathway	2.24E-04	16
	hsa04742	Taste transduction	2.54E-04	9
	hsa04720	Long-term potentiation	2.89E-04	8
	hsa04723	Retrograde endocannabinoid signaling	3.97E-04	12
	hsa00220	Arginine biosynthesis	1.78E-03	4
	hsa04726	Serotonergic synapse	2.39E-03	9
	hsa04964	Proximal tubule bicarbonate reclamation	2.54E-03	4
	hsa04971	Gastric acid secretion	2.98E-03	7
	hsa04925	Aldosterone synthesis and secretion	4.91E-03	7
	hsa04010	MAPK signaling pathway	5.79E-03	14
	hsa04924	Renin secretion	6.20E-03	6
	hsa05034	Alcoholism	6.46E-03	11
	hsa04961	Endocrine and other factor-regulated calcium reabsorption	6.75E-03	5

Figure 4.

The network for the top 50 genes in modules. (A) Lightcyan; (B) black; deep color nodes: downregulated.

MS Proteomic Analysis

Compared with control samples, 503 DEPs (262 upregulated and 241 downregulated) were obtained from GBM samples. GO and KEGG pathway results for these proteins are presented in Figure 5.

Figure 5.

The top 5 gene ontology (GO) (A) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (B) pathways for differentially expressed proteins.

For the 503 differentially expressed proteins, the expression values of 495 genes (corresponding to 481 proteins) were detected at the transcriptome level. Of them, 89 genes were also differentially expressed at the proteome level. Furthermore, 25 genes were upregulated and 62 were downregulated at both the transcriptome and proteome levels (Figure 6A). In addition, 2 genes were downregulated at the transcriptome level and upregulated at the proteome level. The 62 coexpressed genes were mainly enriched in the regulation of synaptic vesicle exocytosis, neural nucleus development, and the regulation of neurotransmitter levels (Figure 6B).

Figure 6.

Association analysis of proteome and transcriptome (A) and clueGO enrichment results of co-down-regulated genes/proteins (B).

The results showed that the TCGA database lacked expression information on immunoglobulin heavy chain constant γ and immunoglobulin κ constant, the difference observed for aldehyde dehydrogenase 1 family member A1 was not significant, and the other 84 genes matched with transcriptome sequencing results.

Survival and Prognostic Analysis

A univariate Cox regression analysis was conducted for GBM samples using the 84 genes that were validated by TCGA database, and the P value of prognostic risk for every gene was obtained. The following 8 prognostic risk factors were identified: heparin sulfate proteoglycan 2 (P = 0.007), protein disulfide isomerase family A member 4 (PDIA4) (P = 0.007), annexin A2 (P = 0.011), laminin subunit γ 1 (P = 0.012), FN1 (P = 0.014), reticulon-1 (P = 0.021), chitinase 3-like 1 (CHI3L1) (P = 0.040), and ANXA1 (P = 0.049). A multivariate Cox regression analysis was conducted for the 8 prognostic risk factors. After optimization, 3 genes (CHI3L1, FN1, and PDIA4) were included. Finally, FN1 was found to be a prognostic risk factor for GBM (Table 4). We grouped FN1 based on gene expression values to construct a KM survival curve between patients in high and low expression groups. The logrank test found that there was a significant difference in survival between the high and low expression groups (Figure 7).

Table 4.

Multivariate Cox Regression Analysis for CHI3L1, FN1, and PDIA4.a

	Coefficient	Exponent (coef)	Se (coef)	Z	P value
CHI3L1	1.35e-06	1.00E + 00	7.43e-07	1.82	0.058
FN1	4.83e-06	1.00E + 00	2.36e-06	2.05	0.040
PDIA4	8.06e-05	1.00E + 00	4.28e-05	1.88	0.060

a Likelihood ratio test = 14.9 on 3 df, P = 0.00187.

Figure 7.

Kaplan-Meier curve analysis between high and low expression groups of FN1.

As shown in Figure 8, the expression of FN1 in GBM samples was significantly higher than that of the control samples (P < 0.01).

Figure 8.

The expression of FN1 in glioblastoma multiforme (GBM) samples compared with control samples (“***” represents P < 0.01).

Discussion

In the present study, both DEGs and differentially expressed proteins in human GBM were identified. Upregulated DEGs were mainly enriched in ECM–receptor interaction and PI3K–Akt signaling pathways, whereas downregulated DEGs were mainly enriched in glutamatergic synapses, calcium signaling, and neuroactive ligand–receptor interactions. In total, 25 genes were upregulated and 62 were downregulated at both the transcriptome and proteome levels. The results of TCGA database searching showed that 84 of these 87 genes matched transcriptome sequencing results, and the Cox regression analysis suggested that FN1 was a prognostic risk factor for GBM. Furthermore, experimental verification showed a significantly higher expression of FN1 in GBM samples compared with control samples.

FN1, a high-molecular weight glycoprotein of the extracellular matrix, plays an important role in embryonic development and wound healing. Altered expression, organization, and degradation of FN1 have been related to the pathologies of fibrosis and cancer. FN1 is also involved in glioma cell motility, cohesion, and tumor angiogenesis. Furthermore, this protein plays an important role in the maintenance of integrin b1 fibronectin receptors in glioma cells. Volyanskyy et al have shown that FN1 is highly upregulated in GBM. Han et al have shown that FN1 was associated with Tax-interacting protein 1-regulated angiogenesis in human GBM. Considering these findings together with the results of the present study, it is very likely that FN1 plays a significant role in GBM development and is an important prognostic risk factor.

Furthermore, the present study showed that upregulated DEGs were mainly enriched in ECM–receptor interaction and PI3K–Akt signaling pathways. The ECM–receptor interaction pathway is involved in important features of malignant gliomas and its downregulation has been related to temozolomide-resistant phenotypes. Hu et al identified 6 KEGG pathways related to the progression of glioma, among which the ECM–receptor interaction pathway plays an important role in the occurrence and development of GBM. Some studies have also reported the role of ECM pathway in the development of other cancers, such as kidney cancer and oral squamous cell carcinoma. On the other hand, PI3K–Akt signaling has also been shown to be involved in malignant transformation, growth, proliferation, and metastasis of ovarian cancer. Li et al have indicated that PI3K–Akt signaling plays an important role in glioma development, but therapeutic targeting of this pathway did not provide satisfactory results. In contrast, other studies have shown that the PI3K–Akt signaling pathway is a therapeutic target for estrogen receptor-positive breast, ovarian, and non-small-cell lung cancers. Therefore, both ECM–receptor interaction and PI3K–Akt signaling pathways are involved in GBM development. Furthermore, results from the current study showed that FN1 was mainly enriched in the ECM–receptor interaction and PI3K–Akt signaling pathways. Considering all the above, we infer that FN1 plays an important role in GBM development via the ECM–receptor interaction and PI3K–Akt signaling pathways.

The results of the present study suggest that FN1 is a significant prognostic risk factor for GBM and may play a role in GBM development and progression through the ECM–receptor interaction and PI3K–Akt signaling pathways. Furthermore, FN1 may be considered as a therapeutic target for GBM in future clinical trials. However, the sample size of this study was small; further studies with larger sample sizes are needed.