Scouting for common genes in the heterogenous hypoxic tumor microenvironment and their validation in glioblastoma

Abstract

Investigating the therapeutic and prognostic potential of genes in the heterogeneous hypoxic niche of glioblastoma. We have analyzed RNA expression of U87MG cells cultured in hypoxia compared to normoxia. Common differentially expressed genes (DEGs) from GSE45301 and GSE18494 and their functional enrichment was performed using MetaScape and PANTHER. Hub genes and their ontology were identified using MCode cytoHubba and ClueGO and validated with GlioVis, Oncomine, HPA and PrognoScan. Using the GEO2R analysis of GSE45301 and GSE18494 datasets, we have found a total of 246 common DEGs (180 upregulated and 66 downregulated) and identified 2 significant modules involved in ribosome biogenesis and TNF signaling. Meta-analysis of key genes of each module in cytoHubba identified 17 hub genes (ATF3, BYSL, DUSP1, EGFR, JUN, ETS1, LYAR, NIP7, NOLC1, NOP2, NOP56, PNO1, RRS1, TNFAIP3, TNFRSF1B, UTP15, VEGFA). Of the 17 hub genes, ATF3, BYSL, EGFR, JUN, NIP7, NOLC1, PNO1, RRS1, TNFAIP3 and VEGFA were identified as hypoxia signatures associated with poor prognosis in Glioma. Ribosome biogenesis emerged as a vital contender of possible therapeutic potential with BYSL, NIP7, NOLC1, PNO1 and RRS1 showing prognostic value.

Supplementary information: The online version contains supplementary material available at 10.1007/s13205-021-02987-2.

Keywords: GEO2R; Glioma; Hub genes; Hypoxia; PPI network; Ribosome biogenesis.

Fig. 1

Volcano plot and Venn diagram of differentially expressed genes. Differentially expressed genes were selected by volcano plot filtering (fold change ≥ 1 and p-value ≤ 0.05) from datasets A GSE18494 and B GSE4530. The Blue dots represents upregulated genes while the red represents downregulated genes. Common differentially C upregulated genes and D downregulated between GSE118494 and GSE45301 were determined. The Blue circle represents upregulated genes while peach represents downregulated genes. A darker tone represents an intersection between the two

Fig. 2

Functional enrichment analysis of the overlapping differentially expressed genes with Metascape. A Bar graph of enriched terms of the upregulated overlapping genes (colored by p values). B The network of enriched terms (colored by cluster ID) of the upregulated genes. C The network of enriched terms (colored by − log 10 p-value) of the upregulated genes. D Bar graph of enriched terms of the down-regulated overlapping genes (colored by p values). E The network of enriched terms (colored by cluster ID) of the down-regulated genes. F The network of enriched terms (colored by − log 10 p-value) of the down-regulated genes. Nodes represent enriched terms or pathways with node size indicating the number of DEGs involved in. Nodes represent enriched terms or pathways or Hallmark gene sets with node size indicating the number of DEGs involved in. Nodes sharing the same cluster are typically close to each other, and the thicker the edge displayed, the higher the similarity

Fig. 3

Overlapping differentially expressed genes associated with functional gene ontology (GO) classification system using PANTHER. Upregulated A biological process B cellular component C molecular function and downregulated D biological process E cellular component F molecular function are represented

Fig. 4

Protein–protein interaction network based on 246 overlapping differentially expressed genes (DEGs) visualized by Cytoscape. The red nodes represent upregulated genes and the green nodes represent downregulated genes. The lines represent the interaction between proteins. Disconnected nodes were hidden in the network

Fig. 5

Significant modules identified from the protein–protein interaction network by Cytoscape A Module1: score = 8.667 and B Module2: score = 5.875. Gene annotation functional enrichment analysis for each module by ClueGo is represented for module 1 as B and C; for module 2 as E and F)

Fig. 6

Employing an online resource, we used five intersecting algorithms to generate a Venn plot to identify significant hub genes from A module 1 and B Module 2. Areas with different colors correspond to different algorithms. The cross areas indicate the commonly accumulated in DEG modules. The elements in concurrent areas are the hub genes. 9 genes were identified in module 1 (PNO1, NOLC1, NIP7, UTP15, BYSL, RRS1, NOP2, NOP56, LYAR) and eight genes in module 2 (VEGFA, JUN, TNFRSF1B, ETS1, ATF3, TNFAIP3, EGFR, DUSP1)

Fig. 7

Interaction network and functional enrichment of the identified hub genes. The nodes indicate the top hub genes, and the edges indicate close interactions between the core genes. The interactions were identified as highly significant as represented with the colour code (darker represents highly significant interactions)

Fig. 8

GSEA enrichment between Normoxia and Hypoxia risk groups in GSE45301 and GSE18494 A Significant Gene Ontology (GO) enrichment plots of Hub genes and B Significant Hallmark enrichment plots of Hub genes. Normalized enrichment score (NES) > 1, nominal p-value (NOM p-value) < 0.05 and False discovery rate (FDR) < 0.25 were considered significant gene sets

Fig. 9

Expression of hub genes in glioblastoma vs normal brain tissue. Statistical significance was indicated in the figures, ***p < 0.001; **p < 0.01; *p < 0.05; ns, not significant. ATF3, BYSL, DUSP1, EGFR, JUN, NIP7, NOP2, NOP56, PNO1, TNFAIP3, TNFRSF1B and VEGFA were significantly upregulated while NOLC1 is significantly downregulated in TCGA GBM dataset

Fig. 10

Expression of the hub genes in high-grade glioma (HGG, Glioblastoma) compared to three subtypes of low-grade gliomas (LGG -Oligodendroglioma, Oligoastrocytoma, Astrocytoma). Statistical significance was indicated in the figures, ***p < 0.001; **p < 0.01; *p < 0.05; ns, not significant. ATF3, BYSL, DUSP1, EGFR, ETS1, JUN, NIP7, NOP2, NOP56, PNO1, RRS1, TNFAIP3, TNFRSF1B and VEGFA were significantly upregulated while NOLC1 was significantly downregulated in TCGA GBM-LGG datasets

Fig. 11

Meta-analysis of differential mRNA expression levels of hub genes in glioblastoma datasets. Seven datasets of glioblastoma samples and normal brain samples in the Oncomine were utilized to determine that ATF3, BYSL, EGFR, ETS1, JUN, NIP7, NOP2, NOP56, PNO1, TNFAIP3, TNFRSF1B and VEGFA was significantly upregulated and increased expression displayed by DUSP1, LYAR, RRS1 and UTP15 was not significant. NOLC1 was observed to be significantly downregulated

Fig. 12

Immunohistochemical analysis of the hub genes. The protein expression of identified hub genes was assessed in glioma cancer tissues and normal glial tissue derived from Human Protein Atlas database. ATF3, BYSL, EGFR, ETS1, JUN, NIP7, RRS1, TNFAIP3, TNFRSF1B and VEGFA) were enhanced in glioma patients

Fig. 13

Association of the expression of hub genes with overall survival prognosis in glioma patients. A Kaplan–Meier curves shows that low expression of BYSL, NOLC1, NIP7, PNO1 and RRS1 as well as higher expression of BYSL and NOLC1 with significant cox p values < 0.05 and hazard ratio (95% CI) of < 1 associated with better prognosis. Higher expression of ATF3, EGFR, JUN, TNFAIP3, and VEGFA with worse prognosis due to hazard ratio (95% CI) of > 1