Identification of differentially expressed key genes between glioblastoma and low-grade glioma by bioinformatics analysis

Data of microarray

Gene expression profile of GSE52009, GSE53733, and GSE4290 were downloaded from GEO database, which is a public and freely accessible database. Based on Agilent GPL6480 platform (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F), GSE52009 dataset included 120 samples, containing 60 WHO II samples and 24 WHO IV samples. GSE53733 was based on the GPL570 platform ((HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array), which contained 70 GBM samples. GSE4290 was based on the GPL570 platform ((HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array), which included 180 samples, containing 76 WHO II samples and 81 WHO IV samples.

Screen genes of differential expression

Differentially expressed genes between low-grade glioma and high-grade glioma was detected by GEO2R, which was an online analysis tool based on R language (Davis & Meltzer, 2007). We set the adjust P-value < 0.05 and |logFC| ≥ 2 as the selection criteria to decrease the false positive rate and false discovery rate. Furthermore, the top 15 genes with higher degree of connectivity were selected as hub genes among the 133 discovered DEGs which includes 56 downregulated genes and 77 upregulated genes. In addition, we used visual hierarchical cluster analysis to show the two groups by Morpheus online analysis software (https://software.broadinstitute.org/morpheus/) and volcano plot of two groups by ImageGP (http://www.ehbio.com/ImageGP/index.php/Home/Index/index.html) after the relative raw data of TXT files were downloaded.

Gene ontology and KEGG pathway analysis of DEGs

With functions including MF, biological pathways, and CC, gene ontology (GO) analysis we annotated genes and gene products (Gene Ontology Consortium, 2006). KEGG comprises a set of genome and enzymatic approaches and abiological chemical energy online database (Kanehisa & Goto, 2000). It is a resource for systematic analysis of gene function and related high-level genome functional information. Database for annotation, visualization and integrated discovery (DAVID) (https://david.ncifcrf.gov/) can provide systematic and comprehensive biological function annotation information for high-throughput gene expression (Dennis et al., 2003). Therefore, we applied GO and KEGG pathway analyses to the DEGs by using DAVID online tools at functional level. We considered P < 0.05 had significant differences. In addition, we used visual analysis to show GO Enrichment plot of two groups by ImageGP (http://www.ehbio.com/ImageGP/index.php/Home/Index/index.html) after the relative raw data of TXT files were downloaded (Geng et al., 2018).

PPI network and module analysis

The online tool, search tool for the retrieval of interacting genes (STRING), is designed to demonstrate the interaction between different proteins (Szklarczyk et al., 2015). STRING in Cytoscape was applied and mapped the DEGs into STRING to detect the possible relationship among the selected DEGs. We set the confidence score ≥0.4, maximum number of interactors = 0 as the selection criteria. In addition, the molecular complex detection (MCODE) was used to screen modules of PPI network in Cytoscape with degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and max. depth = 100. DAVID was used to perform the signal pathway analysis of genes in the module. A total of 15 hub genes were also mapped into STRING with confidence score ≥0.4, maximum number of interactors ≤5. The potential information was explored through GO and KEGG pathway analysis.

Comparing the expression level of the hub genes

GlioVis (http://gliovis.bioinfo.cnio.es/) is a user-friendly web application for data visualization and analysis to explore brain tumors expression datasets, which was used to analyze the gene expression data of brain tumors and normal samples based on the TCGA datasets. (Bowman et al., 2017). The customizable functions are provided such as analyzing the differences of expression levels between glioblastoma and low-grade glioma, so the expression of these genes was demonstrated. And the relationship could be visualized through the boxplot. All values are presented as the mean ± SD. All statistical analyses were performed by SPSS 19.0 software. A difference of P < 0.05 was considered statistically significant.

Gene expression profile and gene set enrichment analysis

The expression profiles of GSE53733 were downloaded from the GEO data base. We used GSEA (http://www.broadinstitute.org/gsea) to detect the potential genes influenced by Serum amyloid A1 (SAA1) and tissue inhibitor of metalloproteinases-1 (TIMP1) through Java programming. According to their hub genes expression level (top 50%: high vs. bottom 50%: low), we divided the patients into two groups, and GSEA was conducted to analyze the effects of selected genes expression level on different BP. We set P-value of <0.05 and false discovery rates of <0.25 as selection criteria to confirm significant gene sets.

Human tissue samples

Low-grade gliomas and GBM tissues were collected from the Department of Neurosurgery, Renmin Hospital of Wuhan University, Wuhan, China. The clinical glioma specimens were examined and diagnosed by pathologists at Renmin Hospital of Wuhan University. This study was approved by the Institutional Ethics Committee of the Faculty of Medicine at Renmin Hospital of Wuhan University (approval number: 2012LKSZ (010) H). Informed consent was obtained from all patients whose tissues were used.

RNA extraction and quantitative real-time PCR

Total RNA from cancer tissues was prepared using Trizol reagent (Invitrogen, Carlsbad, CA, USA), and cDNA was synthesized using a PrimeScript RT Reagent Kit with gDNA Eraser (RR047A; Takara, Kusatsu, Japan). Quantitative real-time PCR (qPCR) for SAA1 and TIMP1 mRNA levels were performed using SYBR Premix Ex Taq II (RR820A, Takara, Kusatsu, Japan) according to the manufacturer’s instructions and performed in Bio-Rad CFX Manager 2.1 real-time PCR Systems (Bio-Rad, Hercules, CA, USA). GAPDH was used as internal controls. The data were analyzed by the relative Ct method and expressed as a fold change compared with the control. The primer sequences included the following: GAPDH 5′-GGAGCGAGATCCCTCCAAAAT-3′(Forward), 5′-GGCTGTTGTCATACTTCTCATGG-3′(Reverse); SAA1 5′-CCTGGGCTGCAGAAGTGATCAGCGA-3′(Forward), 5′-AGTCCTCCGCACCATGGCCAAAGAA-3′(Reverse); TIMP1 5′-CTTCTGCAATTCCGACCTCGT-3′(Forward), 5′-ACGCTGGTATAAGGTGGTCTG-3′(Reverse).

Identification of DEGs and hub genes

A total of 60 WHO II samples and 24 WHO IV samples from GSE52009 were selected in this study. DEGs were detected by applying the GEO2R online analysis tool, setting adjust P-value < 0.05 and |logFC| ≥ 2 as selection criteria. A total of 133 differential expressed genes, containing 77 upregulated genes and 56 downregulated genes, were detected after the analysis of GSE52009. In addition, we selected 15 hub genes with higher degree of connectivity (Table 1). The results were validated with a DEG expression heatmap and volcano plot of the all downregulated genes and upregulated genes (Fig. 1).

GO function and KEGG pathway enrichment analysis

To explore more particular knowledge of the selected DEGs, we used DAVID to gain the results of GO function and KEGG pathway enrichment analysis. All DEGs were imported to DAVID software, and GO analysis results demonstrated that upregulated and downregulated DEGs were particularly enriched in the following biological processes (BP): cell migration, locomotion and leukocyte migration, cell motility for upregulated DEGs, and for downregulated DEGs nervous system development, brain development and regulation of cell projection organization (Table 2; Figs. 2A and 2B). The upregulated DEGs were enriched in phospholipase A2 inhibitor activity, growth factor binding, extracellular matrix (ECM) structural construction, receptor binding, and the downregulated DEGs were enriched in calcium ion binding, structural construction of myelin sheath, and protein complex binding for MF (Table 2; Fig. 2). Moreover, GO CC analysis showed that the upregulated DEGs were enriched in the proteinaceous ECM, ECM and cytoplasmic membrane-bounded vesicle lumen, and downregulated DEGs enriched in neuron part, myelin sheath, and internode region of axon (Table 2; Figs. 2A and 2B).

Interestingly, the most significantly enriched of KEGG pathway only showed in upregulated pathway, including adhesion, ECM-receptor interaction, amoebiasis, and PI3K-Akt signaling pathway (Table 3; Fig. 2C).

Hub genes and module screening from PPI network

Protein–protein interaction network of the top 15 hub genes with higher degree of connectivity was made based on the information in the STRING protein query from public databases (Fig. 3A). The top module was selected by using MCODE plug-in in the PPI network (Fig. 3B).

The Kaplan–Meier plotter of hub genes

The website, http://gepia.cancer-pku.cn/, freely provides the prognostic data of the hub genes. It was found that expression of VEGFA (HR 4.2, P < 0.001) was associated with worse OS for glioblastoma patients, as well as NDC80 (HR 5.8, P < 0.001), CENPA (HR 5.3, P < 0.001), CENPF (HR 3.9, P < 0.001), Non-SMC condensin I complex subunit G (NCAPG) (HR 5.6, P < 0.001), ASPM (HR 5, P < 0.001), ITGA2 (HR 3, P < 0.001), TIMP1 (HR 7, P < 0.001)and SAA1 (HR 4.8, P < 0.001) (Fig. 4).

Expression level and relationship with molecular pathologic diagnosis of hub genes

We used data from GlioVis to detect the hub gene expression level between GBM and LGG including astrocytoma, oligodendroglioma, and oligoastrocytoma, the expression level of SAA1 and TIMP1 significantly increased in GBM (Figs. 5A and 5C). The expression levels of SAA1 have no significant difference in three kind of LGG (Fig. 5B). However, the expression level of TIMP1 is significantly higher in astrocytoma than oligodendroglioma and oligoastrocytoma (Fig. 5D). We further verified our finding in the GSE4290 dataset and got consistent result (Fig. S1). Then we detected the sample collected in our hospital and found both SAA1 and TIMP1 are significantly increased in GBMs compared with LGG (Fig. S2; Table S1). We also detect the relationship between expression level and molecular pathologic diagnosis of hub genes. We found both SAA1 and TIMP1 increase in both Isocitrate dehydrogenase (IDH) mutant IDH wild type. The same results were found in MGMT promoter and non-deletion of chromosome 1p.19q. Because of the limited samples in the datasets, we didn’t the result of co-deletion of chromosome 1p.19q. (Figs. 5E and 5F) Further, we also found that both SAA1 and TIMP1 played important roles in MES-like in the Isocitrate dehydrogenase (IDH) wild type (Figs. 5G and 5H).

Gene expression profile and gene set enrichment analysis

We managed a GSEA by using GBM patient gene profiling data (GSE53733), and showed in Fig. 5, gene set differences in SAA1 in low versus high glioma patients indicated that SAA1 regulates biology process mainly associated with inflammatory response processes (P < 0.001 FDR = 0.012) and cytokine mediated signaling pathway (P < 0.001 FDR = 0.012) (Figs. 5I and 5J). We considered that SAA1 may negatively regulate inflammatory response and might promote the survival of cancer cells. We concluded that TIMP1 might negatively regulates adaptive immune response based on somatic recombination of immune receptors built from a leucine-rich superfamily (P < 0.001 FDR = 0.021) and response to interferon (P < 0.001 FDR = 0.027) may promote the survival of cancer cells (Figs. 5K and 5L).