The present study aimed to identify differentially expressed genes (DEGs) regulated by transcription factors (TFs) in glioblastoma, by conducting a bioinformatics analysis. The results of the present study may provide potential therapeutic targets that are involved in the development of glioblastoma. The GSE4290 raw data set was downloaded from the Gene Expression Omnibus database, and consisted of 23 non-tumor samples and 77 glioblastoma (grade 4) tumor samples. Robust Multichip Averaging was used to identify DEGs between the glioblastoma and non-tumor samples. Functional enrichment analysis of the DEGs was also performed. Based on the TRANSFAC® database, TFs associated with the glioblastoma gene expression profile were used to construct a regulatory network. Furthermore, trimmed subnets were identified according to calculated Z-scores. A total of 676 DEGs were identified, of which 190 were upregulated and 496 were downregulated. Gene Ontology analysis demonstrated that the majority of these DEGs were functionally enriched in synaptic transmission, regulation of vesicle-mediated transport and ion-gated channel activity. In addition, the enriched Kyoto Encyclopedia of Genes and Genomes pathway included neuroactive ligand-receptor interaction, calcium signaling pathway, p53 signaling pathway and cell cycle. Based on the TRANSFAC® database, transcriptional regulatory networks with 2,246 nodes and 4,515 regulatory pairs were constructed. According to the Z-scores, the following candidate TFs were identified: TP53, SP1, JUN, STAT3 and SPI1; alongside their downstream DEGs. TP53 was the only differentially expressed TF. These candidate TFs and their downstream DEGs may have important roles in the progression of glioblastoma, and could be potential biomarkers for clinical treatment.
Glioblastoma is the most frequent and aggressive brain malignancy in adults, and is characterized by a heterogeneous population of cells that are involved with progression of the disease (1). It is a rapidly fatal malignancy and the majority of patients with glioblastoma suffer from a poor quality of life (2,3). Currently, the standard clinical treatment is surgical resection of the malignant tissues, followed by radiotherapy and chemotherapy (4–7). However, patients that receive these treatments may rapidly develop resistance to chemotherapy (8). Recent studies have focused on the identification of candidate biomarkers of glioblastoma development, in order to produce a more effective therapeutic strategy (9–11).
Transcription factors (TFs) have important roles in the transcriptional networks that regulate gene expression, and modify and control cancer phenotypes (12,13). Differentially expressed TFs in glioblastoma, and their downstream gene targets, may be potential therapeutic biomarkers of glioblastoma (12,13). O6-methylguanine DNA methyltransferase (MGMT) promoter hypermethylation (14,15) and isocitrate dehydrogenase 1 (16–18) have previously been suggested as potential therapeutic targets, and regulation of MGMT expression has been reported in numerous clinical studies (19,20). It has been suggested that MGMT expression may be regulated by inhibiting its upstream TF, such as SP1 in glioblastoma (21).
Sun et al (22) collected mRNA expression data (GSE4290) from patients with brain tumors, and demonstrated that downregulation of stem cell factor (SCF) inhibits tumor-mediated angiogenesis and glioma growth in vivo, whereas overexpression of SCF was associated with reduced survival in patients with malignant glioma. Numerous studies have identified glioblastoma-associated genes based on the GSE4290 dataset, with the aim of improving diagnosis of glioma at the molecular level (23,24). However, the importance of differentially expressed TFs has yet to be explored. The present study aimed to identify the differentially expressed TFs in glioblastoma, and the corresponding critical pathways involved in glioblastoma development.
In the present study, the raw mRNA data of Sun et al (22) was downloaded from the Gene Expression Omnibus (GEO), and the differentially expressed genes (DEGs) in glioblastoma samples were identified. Functional enrichment analysis of the DEGs was then performed. TFs associated with the glioblastoma gene expression profile were used to construct a regulatory network. The present study may improve understanding regarding the development of glioblastomas. Furthermore, the differentially expressed TFs may be potential biomarkers for the prognosis and therapy of glioblastoma.
Databases and methodsData acquisition
The raw data was downloaded from the GSE4290 dataset (22) deposited in the GEO (http://www.ncbi.nlm.nih.gov/geo/)(25). The dataset included 23 samples from patients with epilepsy, which are considered non-tumor samples, and 77 glioblastoma (grade 4) tumor samples. The platform was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array.
Analysis of DEGs
The raw data was initially analyzed using R software (v.3.0.0; http://www.r-project.org/). The chip data was normalized using the Robust Multichip Averaging method (26) in Affy package (http://www.r-project.org/) (27). The DEGs were then identified using the Limma package (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (28) and tested for multi-test correction by Bayes law (29). Genes with P<0.05 and |log2fold change (FC)| >1.5 were considered to be DEGs between the tumor and non-tumor groups.
Functional enrichment analysis
For functional analysis of the selected DEGs, the DEGs were imported into the Database for Annotation, Visualization and Integrated Discovery (http://david.abcc.ncifcrf.gov/) (30), in order to perform a Gene Ontology (GO) functional enrichment analysis and a Kyoto Encyclopedia of Genes and Genomes (KEGG) (31,32) pathway enrichment analysis. GO analysis encompasses three domains: Biological process, cellular components and molecular functions. P<0.05 was considered to indicate significance.
Weight of regulatory network
Based on the TRANSFAC® (33) database (http://www.gene-regulation.com/pub/databases.html) and the glioblastoma gene expression profile (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE4290), TFs identified in the two datasets were selected and used to establish a regulatory network with their target genes. Combined with the gene expression levels, formulae i and ii were used to calculate the average rank correlation coefficient and formula iii was used to calculate the difference value of Spearman coefficients within the regulatory network. The absolute values of the averages of rank correlation coefficient were defined as weight of TF-gene pairs and the absolute value of difference value was defined as weighted coefficient (28).
where Eij is the TF-target gene between TF Vi and gene Vj; k is the kth sample; Vi and Vj are ranked by their expression levels in the samples respectively, and Xjk is the rank of Vi in kth sample, Xik is the rank of Vj of kth sample; xi, xj are the average ranks of Vi and Vj in the samples, respectively. rEij1 and rEij2 represent the Spearman coefficients of Eij in compared samples respectively. Permutation test was applied to rank the random difference values. TF-gene pairs with a weighted coefficient >90% of the weighted coefficient value were excluded from further analysis (34).
Screening of sub-networks within the regulatory network
TFs with a degree >15 in the regulatory network were selected and used to establish sub-networks with their target genes. The weight of TF-gene pairs in the sub-networks were scored using the following methods. Initially, the weighted coefficients of all TF-gene pairs within the regulatory network were ranked and defined as a background set (E), whereas the sub-networks were considered as an objective set (S). The score of S enriched into E was then calculated by gene set enrichment analysis (35), according to formula iv:
where Ej is the jth TF-target in the ranked regulatory pairs; rj is the weight of the jth regulatory pair in background set; P is a parameter and set as 1; N is the number of regulatory pairs in E; NH is the number of regulatory pairs in the subnet S. The enrichment score (ES) is the maximum deviation between Phit and Pmiss.
TF-gene pairs without contribution to the ES were excluded from the analysis (34,35). To estimate the significance of ES of the sub-regulatory networks, ES was converted into Z value (34) using formula v.
Zs=ES-ES¯S′
where ES (bar) is the mean of the random ES set; and S’ is the standard deviation of the random ES set.
DEGs in the trimmed subnet
Genes from the gene expression profile were defined as a background set (E), whereas genes in the trimmed subnet were defined as an objective set (S). The P-values of the DEGs were ranked and the ES was calculated using formula vi (34).
where gj is the jth gene in the ranked genes; rj is the magnitude of differential expression of the jth gene; P is a parameter and set as 1; M is the number of genes in L; and MH is the number of genes in Strimmed.
Genes that did not contribute to the trimmed subnet ES were excluded. The significance of DEGs in the trimmed subnet was calculated by Z value transformation. The final Z-score was calculated using formula vii. The top five trimmed subnets were selected as the candidate regulatory subnets in glioblastoma.
ZScombinedi=Zsi′+Ztrimmedi′ResultsIdentification of DEGs and functional enrichment analysis
With a cut-off value of P<0.05 and |log2FC| >1.5, a total of 676 DEGs were identified, of which 190 were upregulated and 496 were downregulated (Table I). GO analysis demonstrated that the majority of DEGs were enriched in synaptic transmission, regulation of vesicle-mediated transport and ion-gated channel activity (Fig. 1). In addition, KEGG pathway enrichment analysis identified the significantly enriched pathways, which included neuroactive ligand-receptor interaction, calcium signaling pathway, p53 signaling pathway and cell cycle (Fig. 2).
Establishment of a weighted regulatory network and trimmed subnets
To identify TFs in the DEGs, TF-gene pairs were selected based on the TRANSFAC® database and a transcriptional regulatory network (not weighted) with 2,246 nodes and 4,515 regulatory pairs was constructed (Fig. 3). With a weighted coefficient >90% of the random weighted coefficient, 1,312 pairs were excluded by permutation test.
TF-gene pairs of trimmed subnets were calculated and the corresponding DEGs were scored. According to the Z-scores, genes with the top 10 highest Z-scores were identified and the corresponding subnets were constructed (Fig. 4). The candidate TFs and their downstream DEGs are listed in Table II. Only TP53 was identified as a differentially expressed TF in glioblastoma.
Discussion
In order to identify potential biomarkers for glioblastoma prognosis and therapy, a bioinformatics analysis was performed on the GSE4290 dataset. A total of 676 DEGs were identified, of which 190 were upregulated and 496 were downregulated. The majority of DEGs were functionally enriched in synaptic transmission, regulation of vesicle-mediated transport and ion-gated channel activity. Furthermore, the enriched KEGG pathways of DEGs included neuroactive ligand-receptor interaction, calcium signaling pathway, p53 signaling pathway and cell cycle. Based on the TRANSFAC® database, a transcriptional regulatory network consisting of 2,246 nodes and 4,515 regulatory pairs was constructed. Based on weighted Z-scores, TP53, SP1, JUN, STAT3, and SPI1 were identified as crucial TFs involved in the development of glioblastoma.
As a common malignancy with poor prognosis, glioblastoma tumors harbor various cell types, including vascular cells, microglia, peripheral immune cells and neural precursor cells, which indicates that there is active communication ongoing between the tumor cells and non-tumor cells, and there is a dramatic turnover in the microenvironment (1). It has previously been shown that calcium-mediated transduction systems, together with active gap junctions, have key roles in the communication of GL15 human glioblastoma cells with surrounding cells (36). Eukaryotic cells are capable of using multivesicular bodies for cytoplasmic trafficking and release of exosomes, which may transfer genetic information between non-immune cells (37). The pathway enrichment results of the present study demonstrated that the DEGs in glioblastoma were enriched in ion-gated channels, gap junction signaling, vesicle-mediated transport signaling and cell-cell signaling, supporting the crucial role of calcium transport, coupled with gap junctions, in the invasive capabilities of glioblastoma. Therefore monitoring these signaling pathways may aid prediction of tumor progression.
Among the five TFs identified in the present study to be associated with glioblastoma, TP53 was the only DEG. TP53 encodes p53, a well-known tumor suppressor protein (38). The abnormal expression of p53 leads to failures of cell cycle and apoptosis regulation, as well as cancer development (38). However, few studies have investigated the role of p53 as a TF. Notably, the present study also identified JUN as a candidate biomarker, which is a proto-oncogene that encodes a component of the mitogen-inducible immediate-early TF AP1 and c-Jun, and regulates the cell cycle (39). It has previously been reported that the regulation of JUN in the cell cycle and apoptosis is associated with p53 (40). Furthermore, overexpression of MGMT has previously been shown to accompany an increased recruitment of c-Jun in glioblastoma (20); however, the association between TP53 and JUN in glioblastoma progression has yet to be elucidated. TP53 and JUN may act as potential biomarkers for the prognosis of glioblastoma.
SP1 was also identified as a candidate TF and the majority of its downstream targets were differentially expressed in glioblastoma, thus indicating that SP1 may be critical for the development of glioblastoma. Previous studies (19,20) have targeted the transcriptional activity of SP1 to regulate the expression of MGMT and other genes for glioblastoma therapy. SPI1 is also a putative proto-oncogene associated with tumor progression (41), which encodes a protein that functions in the development of lymphocytes (42). SPI1 may influence the development of glioblastoma through regulation of functional immune cells. SPI1 may also be a potential biomarker or therapeutic target for glioblastoma; however, this requires further confirmatory study.
In conclusion, the present study identified DEGs between glioblastoma and non-tumor samples, and a functional enrichment analysis of the DEGs was performed. According to Z-scores, the candidate TFs: TP53, SP1, JUN, STAT3 and SPI1, and their downstream DEGs, may have important roles in the progression of glioblastoma, and may be potential biomarkers for clinical treatment.
ReferencesCharlesNAHollandECGilbertsonRGlassRKettenmannHThe brain tumor microenvironmentGlia60502514201210.1002/glia.2126422379614BucknerJCFactors influencing survival in high-grade gliomasSemin Oncol301014200310.1053/j.seminoncol.2003.11.031OmuroADeAngelisLMGlioblastoma and other malignant gliomas: a clinical reviewJAMA31018421850201310.1001/jama.2013.28031924193082LinzUCommentary on Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial (Lancet Oncol 2009; 10: 459–466)Cancer11618441846201010.1002/cncr.2495020151424StuppRHegiMEMasonWPEffects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trialLancet Oncol10459466200910.1016/S1470-2045(09)70025-719269895MinnitiGScaringiCBaldoniAHealth-related quality of life in elderly patients with newly diagnosed glioblastoma treated with short-course radiation therapy plus concomitant and adjuvant temozolomideInt J Radiat Oncol Biol Phys86285291201310.1016/j.ijrobp.2013.02.01323642624OikeTSuzukiYSugawaraKRadiotherapy plus concomitant adjuvant temozolomide for glioblastoma: Japanese mono-institutional resultsPLoS One8e78943201310.1371/journal.pone.0078943242657313827088WenPYKesariSMalignant gliomas in adultsN Engl J Med359492507200810.1056/NEJMra070812618669428Mehrian-ShaiRChenCShiTInsulin growth factor-binding protein 2 is a candidate biomarker for PTEN status and PI3K/Akt pathway activation in glioblastoma and prostate cancerProc Natl Acad Sci USA10455635568200710.1073/pnas.0609139104173722101838515SorensenAGBatchelorTTZhangWTA “vascular normalization index” as potential mechanistic biomarker to predict survival after a single dose of cediranib in recurrent glioblastoma patientsCancer Res6952965300200910.1158/0008-5472.CAN-09-0814195498892824172TitulaerMKMustafaDASiccamaIA software application for comparing large numbers of high resolution MALDI-FTICR MS spectra demonstrated by searching candidate biomarkers for glioma blood vessel formationBMC Bioinformatics9133200810.1186/1471-2105-9-133183126842323386PapagiannakopoulosTShapiroAKosikKSMicroRNA-21 targets a network of key tumor-suppressive pathways in glioblastoma cellsCancer Res6881648172200810.1158/0008-5472.CAN-08-130518829576ChengLBaoSRichJNPotential therapeutic implications of cancer stem cells in glioblastomaBiochemical Pharmacol80654665201010.1016/j.bcp.2010.04.035HamiltonMGRoldánGMaglioccoADetermination of the methylation status of MGMT in different regions within glioblastoma multiformeJ Neurooncol102255260201110.1007/s11060-010-0307-5WangYChenXZhangZComparison of the clinical efficacy of temozolomide (TMZ) versus nimustine (ACNU)-based chemotherapy in newly diagnosed glioblastomaNeurosurg Rev377378201410.1007/s10143-013-0490-xLewandowskaMAFurtakJSzylbergTAn analysis of the prognostic value of IDH1 (isocitrate dehydrogenase 1) mutation in Polish glioma patientsMol Diagn Ther184553201410.1007/s40291-013-0050-73899509SansonMMarieYParisSIsocitrate dehydrogenase 1 codon 132 mutation is an important prognostic biomarker in gliomasJ Clin Oncol2741504154200910.1200/JCO.2009.21.983219636000YanHParsonsDWJinGIDH1 and IDH2 mutations in gliomasN Engl J Med360765773200910.1056/NEJMoa0808710192286192820383SeznecJSilkenstedtBNaumannUTherapeutic effects of the Sp1 inhibitor mithramycin A in glioblastomaJ Neurooncol101365377201110.1007/s11060-010-0266-xKitangeGJMladekACCarlsonBLInhibition of histone deacetylation potentiates the evolution of acquired temozolomide resistance linked to MGMT upregulation in glioblastoma xenograftsClin Cancer Res1840704079201210.1158/1078-0432.CCR-12-0560226751723716364BocangelDSenguptaSMitraSBhakatKKp53-Mediated down-regulation of the human DNA repair gene O6-methylguanine-DNA methyltransferase (MGMT) via interaction with Sp1 transcription factorAnticancer Res29374137502009198469042814523SunLHuiAMSuQNeuronal and glioma-derived stem cell factor induces angiogenesis within the brainCancer Cell9287300200610.1016/j.ccr.2006.03.00316616334KotliarovYKotliarovaSCharongNCorrelation analysis between single-nucleotide polymorphism and expression arrays in gliomas identifies potentially relevant target genesCancer Res6915961603200910.1158/0008-5472.CAN-08-2496191903412644341LiuZXieMYaoZThree meta-analyses define a set of commonly overexpressed genes from microarray datasets on astrocytomasMol Neurobiol47325336201310.1007/s12035-012-8367-5BarrettTWilhiteSELedouxPNCBI GEO: archive for functional genomics data sets - updateNucleic Acids Res41D991D995201310.1093/nar/gks1193IrizarryRAHobbsBCollinFExploration, normalization, and summaries of high density oligonucleotide array probe level dataBiostatistics4249264200310.1093/biostatistics/4.2.24912925520GautierLCopeLBolstadBMIrizarryRAaffy - analysis of Affymetrix GeneChip data at the probe levelBioinformatics20307315200410.1093/bioinformatics/btg40514960456DibounIWernischLOrengoCAKoltzenburgMMicroarray analysis after RNA amplification can detect pronounced differences in gene expression using limmaBMC Genomics7252200610.1186/1471-2164-7-252170296301618401DudoitSGilbertHNvan der LaanMJResampling-based empirical Bayes multiple testing procedures for controlling generalized tail probability and expected value error rates: focus on the false discovery rate and simulation studyBiom J50716744200810.1002/bimj.200710473189321384130579Huang DaWShermanBTTanQThe DAVID Gene Functional Classification Tool: a novel biological module-centric algorithm to functionally analyze large gene listsGenome Biol8R183200710.1186/gb-2007-8-9-r183177849552375021KanehisaMGotoSSatoYFurumichiMTanabeMKEGG for integration and interpretation of large-scale molecular data setsNucleic Acids Res40D109D114201210.1093/nar/gkr9883245020KanehisaMGotoSKEGG: kyoto encyclopedia of genes and genomesNucleic Acids Res282730200010.1093/nar/28.1.27WingenderEThe TRANSFAC project as an example of framework technology that supports the analysis of genomic regulationBrief Bioinform9326332200810.1093/bib/bbn01618436575WuCZhuJZhangXIntegrating gene expression and protein-protein interaction network to prioritize cancer-associated genesBMC Bioinformatics13182201210.1186/1471-2105-13-182228389653464615SubramanianATamayoPMoothaVKGene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profilesProc Natl Acad Sci USA1021554515550200510.1073/pnas.0506580102161995171239896MariggioMAMazzoleniGPietrangeloTCalcium-mediated transductive systems and functionally active gap junctions in astrocyte-like GL15 cellsBMC Physiol14200110.1186/1472-6793-1-41138451032183SkinnerAMO’NeillSLKurrePCellular microvesicle pathways can be targeted to transfer genetic information between non-immune cellsPLoS One4e6219200910.1371/journal.pone.0006219195934432704871KastanMBCanmanCELeonardCJP53, cell cycle control and apoptosis: implications for cancerCancer Metastasis Rev14315199510.1007/BF006902077606818LiebermannDGregoryBHoffmanBAP-1 (Fos/Jun) transcription factors in hematopoietic differentiation and apoptosisInt J Oncol12685138519989472112SchreiberMKolbusAPiuFControl of cell cycle progression by c-Jun is p53 dependentGenes Dev13607619199910.1101/gad.13.5.60710072388316508RayDCulineSTavitainAMoreau-GachelinFThe human homologue of the putative proto-oncogene Spi-1: characterization and expression in tumorsOncogene566366819901693183BusslingerMTranscriptional control of early B cell developmentAnnu Rev Immunol225579200410.1146/annurev.immunol.22.012703.10480715032574
Top 10 GO terms enriched by DEGs. (A) Biological processes; (B) Cellular components; (C) Molecular function. The horizontal axis represents the count of enriched DEGs. The vertical axis represents the different GO terms. GO, gene ontology, DEG, differentially expressed gene; GABA, γ-aminobutyric acid.
The KEGG pathway enrichment analysis of DEGs. The horizontal axis represents the count of enriched DEGs. The vertical axis represents the different KEGG pathways. KEGG, Kyoto Enclyclopedia of Genes and Genomes; DEG, differentialy expressed gene; ECM, extracellular matrix.
Transcriptional regulation network of genes from the gene expression profile of glioblastomas. The triangles represent differentially expressed transcription factors (TFs) (red, differential expression; green, expression without difference). The yellow circles represent differentially expressed genes (DEGs) regulated by TFs (yellow, DEGs; blue, expression without difference).
Corresponding sub-regulation network of the transcriptional regulation network. The triangles represent differentially expressed transcription factors (TFs) (red, differential expression; green, expression without difference). The yellow circles represent differentially expressed genes (DEGs) regulated by TFs (yellow, DEGs; blue, expression without difference).
Top 10 up- and downregulated differently expressed genes (DEGs) in glioblastoma tissue samples.