Mining TCGA database for genes of prognostic value in glioblastoma microenvironment

1Glioblastoma multiforme (GBM) is one of the most fatal brain tumors with a mean survival rate of 35.7% at one Year, 4.7% at five years, and median overall survival (OS) of 14.6 months. To better understand the impacts of genetic composition of tumor on clinical prognosis, comprehensive genome-wide gene expression collections such as The Cancer Genome Atlas (TCGA) have been established to categorize and discover genomic abnormalities in large cohorts across the world. In the TCGA database, according to global gene expression profiles, GBM was initially classified into four subtypes: proneural, neural, classical, and mesenchymal. Of these subtypes, the neural subtype is no longer recognized as a major one due to its lack of tumor-intrinsic characteristics based on several recent reports. In 2016, the updated World Health Organization (WHO) classification integrated molecular parameters with histology and divided GBM into three subtypes: (1) IDH-wildtype, (2) IDH-mutant, and (3) NOS (not otherwise specified). With these progresses, gene expression profiling has been increasingly incorporated with and accepted by clinical diagnostic criteria.

2Tumor cell intrinsic genes especially master transcription factors dictate the initiation, progression, and evolution of GBM. On the other hand, tumor microenvironment has also been reported to critically influence gene expression of tumor tissues, hence the clinical outcomes. Tumor microenvironment is the cellular milieu where the tumor is located. It consists of immune cells, mesenchymal cells, endothelial cells, along with inflammatory mediators and extracellular matrix (ECM) molecules. In the tumor microenvironment, immune and stromal cells are two major types of non-tumor components and have been proposed to be valuable for diagnostic and prognostic assessment of tumors. Algorithms have been developed to predict tumor purity using gene expression data from the TCGA database. For instance, Yoshihara et al designed an algorithm called ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data). In this algorithm, the authors calculated immune and stromal scores to predict the infiltration of non-tumor cells, by analyzing specific gene expression signature of immune and stromal cells. Subsequent reports have soon applied the ESTIMATE algorithm to prostate cancer, breast cancer, and colon cancer, showing the effectiveness of such big-data based algorithms, although utility on immune and/or stromal scores of GBM has not been investigated in detail.

3For the first time in this current work, by taking advantage of both TCGA database of GBM cohorts and ESTIMATE algorithm-derived immune scores, we extracted a list of microenvironment associated genes that predict poor outcomes in GBM patients. Importantly, we have validated such correlation in a different GBM cohort available from the Chinese Glioma Genome Atlas (CGGA) database.

4We downloaded gene expression profiles and clinical information of all 417 GBM patients with initial pathologic diagnosis made between 1989 and 2011 from the TCGA database. Among them, 165 (39.6%) patients were female, 248 (59.5%) cases were male, 4 (0.96%) patients were of unknown gender. Pathological diagnosis included 128 (30.7%) cases of classical subtype, 122 (29.3%) mesenchymal subtype, 64 (15.3%) neural subtype, and 103 (24.7%) cases of proneural subtype. Although the neural subtype was later recommended not to be listed as a major subtype, all GBM cases with complete gene expression data and clinical information in the TCGA were included in our analysis. Based on ESTIMATE algorithm, stromal scores ranged from -3,055.72 to 2,016.62, and immune scores were distributed between -1,448 to 3,210.47, respectively (Figure 1A, B). The average immune scores of mesenchymal subtype cases ranked the highest of all 4 subtypes, followed by that of neural subtype, and classical subtype. The proneural subtype cases had the lowest immune scores (Figure 1A, p < 0.0001). Similarly, the rank order of stromal scores across GBM subtypes from highest to lowest is mesenchymal > neural > classical > proneural (Figure 1B, p < 0.0001), indicating that both immune scores and stromal scores are meaningful in the correlation of subtype classification.

5Based on the 2016 WHO classification, mutation of IDH1 predicts a better prognosis in GBM. We plotted the distribution of immune scores and stromal scores based on the status of IDH1 mutation in GBM cases. IDH1 mutant cases have lower immune scores and stromal scores, although statistically not significant. (Figure 1C, D).

6To find out the potential correlation of overall survival with immune scores and/or stromal scores, we divided the 417 GBM cases into top and bottom halves (high vs. low score groups) based on their scores. Kaplan-Meier survival curves (Figure 1E) showed that median overall survival of cases with the low score group of immune scores is longer than the cases in the high score group (442 d vs. 394 d, p = 0.0537 in log-rank test). Consistently, cases with lower stromal scores also showed longer median overall survival compared to patients with higher stromal scores (Figure 1F, 442 d vs. 422 d, p= 0.1262 in log-rank test), although it was not statistically significant.

7To reveal the correlation of global gene expression profiles with immune scores and/or stromal scores, we compared Affymetrix microarray data of all 417 GBM cases obtained in TCGA database. Heatmaps in Figure 2 showed distinct gene expression profiles of cases belong to high vs. low immune scores/stromal scores groups. For comparison based on immune scores, 480 genes were upregulated and 127 genes downregulated in the high score than the low score group (fold change >1.5, p < 0.05). Similarly, for the high and low groups based on stromal scores, 380 genes were upregulated and 25 genes were downregulated in the high score group (fold change >1.5, p < 0.05). Moreover, Venn diagrams (Figure 2C, D) showed that 374 genes were commonly upregulated in the high scores groups, and 25 genes were commonly downregulated. It is worth mentioning that the differentially expressed genes (DEGs) extracted from the comparison of high vs. low immune scores groups covered the majority of genes extracted from the comparison based on stromal scores. Thus, we decided to focus on these DEGs for all subsequent analysis in this manuscript (Figure 2, Supplementary Figure S1).

8To outline the potential function of the DEGs, we performed functional enrichment analysis of the 480 genes (Supplementary Table S1) upregulated in high-immune scores group. Top gene ontology (GO) terms identified included extracellular matrices, immune and inflammatory response, and chemokine activities and integrin binding (Figure 2E-G).

9To explore the potential roles of individual DEGs in overall survival, we generated Kaplan-Meier survival curves from TCGA database. Among the 480 DEGs upregulated in the high-immune scores group, a total of 258 DEGs (Supplementary Table S2) were shown to significantly predict poor overall survival in log-rank test (p<0.05, selected genes are shown in Figure 3).

10To better understand the interplay among the identified DEGs, we obtained protein-protein interaction (PPI) networks using the STRING tool. The network was made up of 8 modules, which included 224 nodes and 1,282 edges. We selected the top three significant modules for further analysis (Figure 4). For the convenience of description, we named these modules IL6, TIMP1, and TLR2 modules, respectively. In the IL6 module (Figure 4A), 83 edges involving 26 nodes were formed in the network, IL6, IL8, ITGB2, ICAM1, CSF1R, IL1B, and CD163 were the remarkable nodes, as they had the most connections with other members of the module. In the TIMP1 module (Figure 4B), TIMP1, CCR5, CXCL12, SERPINE1, SERPING1, C3AR1, SRGN, and SERPINA3 had higher degree values. For the TLR2 module (Figure 4C), several immune response critical genes occupied the center of the module including TLR2, CCL2, CCL5, IGSF6, and CD14.

11Consistent with PPI network analysis, functional enrichment clustering of these genes showed strong association with immune response as well. A total of 30 GO terms of biological process, 12 GO terms of cellular component, and 5 GO terms of molecular function were identified to be significant (false discovery rate, or FDR<0.05, -log FDR> 1.301). Top GO terms included extracellular exosome and ECM (Figure 5A), immune/inflammatory response, and chemotaxis (Figure 5B), and integrin and proteoglycan binding (Figure 5C). In addition, all the pathways that were yielded from the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (Figure 5D) were associated with immune response.

12To find out whether the genes identified from the TCGA database also are of prognostic significance in additional GBM cases, we downloaded and analyzed gene expression data of a cohort of 114 GBM cases from CGGA, an independent glioma database. A total of 44 genes were validated (Figure 6) to be significantly linked to poor prognosis (Table 1), of which 21 genes were of particular interest as they have not been previously associated with poor outcomes in GBM patients.

13In the current work, we attempted to identify tumor microenvironment related genes that contribute to GBM overall survival in the TCGA database. In particular, by comparing global gene expression in a large number of cases with high vs. low immune scores, we extracted 258 genes involved in extracellular matrix and immune response. Importantly, we were able to validate 44 genes in GBM patients from CGGA, a separate GBM database (Figure 7).

14First, we analyzed 480 differentially expressed genes yielded from comparison of high vs. low immune scores (or stromal scores) groups, and found that many of them were involved in tumor microenvironment, as shown by GO term analysis (Figure 2). This is consistent with previous reports that the functions of immune cells and ECM molecules are interrelated in building tumor microenvironment in GBM.

15Next, we performed overall survival analysis of these 480 genes and identified that 258 were associated with poor outcomes in GBM patients. Moreover, we were able to construct 8 protein-protein interaction modules (Figure 4), all of which were related to immune/inflammation response. Highly interrelated nodes in the modules, including IL6, TIMP1, and TLR2, have been reported to promote proliferation, angiogenesis, migration, and invasiveness in GBM cell lines or patient samples, indicating poor prognosis.

16Finally, by cross validation with CGGA (Figure 6), an independent cohort of 114 GBM patients, we identified 44 tumor microenvironment related genes that showed significant correlation between gene expression and prognosis (Table 1). Of the 44 genes identified, 23 genes (e.g., IL6, IL8, TIMP1, CCR5, SREPINE1, and SREPING1) have been reported to be involved in GBM pathogenesis or significant in predicting overall survival, indicating that our big data-based analysis using TCGA and CGGA cohorts has prognostic values. The remaining 21 genes have not previously been linked with GBM prognosis, and could serve as potential biomarkers for GBM. These include complement encoding genes C3AR1 and C5AR1, TNF superfamily members TNFAIP2 and TNFAIP6, interleukins IL1R2 and IL7R, ECM components COL6A2 and PCOLCE, and cell adhesion molecules ICAM1 and ITGB2.

17We are particularly interested in ITGB2 and ICAM1. From the protein-protein interaction network, ITGB2 and ICAM1 are highly interconnected nodes (Figure 4). ITGB2 encodes integrin subunit β2, which acts as a cell’s mechanical anchor to the ECM. Interestingly, single nucleotide polymorphisms of ITGB2 have been shown to be associated with risk of glioma. As an immunoglobulin supergene family member, ICAM1 is shown to be critically involved in adhesion of cancer cells to ECM, in particular, overexpression of ICAM1 correlates with increased tumor malignancy and poor outcome in lung cancer, clear cell renal cell carcinoma, and mouse GBM model.

18Significant progress has been made on the correlation of overall survival with gene expression in GBMs. Many of these experiments were done in animal tumor models, in vitro tumor cell lines, or small cohorts of patients’ tumor samples. However, the complexity of GBM and GBM microenvironment demands more comprehensive analysis consisting of larger cohorts. Fortunately, the rapid development of whole-genome sequencing has allowed high-throughput tumor databases, including TCGA and CGGA (Chinese Glioma Genome Atlas), to be developed and freely available to the research community. These resources have provided a solid foundation for big data analysis of large GBM cohorts.

19The interplay of GBM and its tumor microenvironment critically affects tumor evolution, which subsequently impacts subtype classification, recurrence, drug resistance, and the overall prognosis of patients. Previous reports have provided elegant analysis on how the activation of tumor-intrinsic genes shapes tumor microenvironment. In the current work, we focused on genes characteristic of microenvironment, which in turn affect the development of GBM and hence contribute to patients’ overall survival. Our results may provide additional data in decoding the complex interaction of tumor and tumor environment in GBM.

20In summary, from functional enrichment analysis of TCGA database applied by ESTIMATE algorithm-based immune scores, we extracted a list of tumor microenvironment related genes. These genes were validated in an independent GBM cohort and that could be useful for outlining the prognosis of GBM patients. Some of the previously ignored genes have the potential to become additional biomarkers for GBM. In addition, it would be extremely interesting to test if this new set of genes, when combined, provide a strong predictor of survival than individual genes. Finally, further investigation of these genes could lead to novel insights into the potential association of tumor microenvironment with GBM prognosis in a comprehensive manner.

21Level 3 gene expression profile (level 3 data) for GBM patients was obtained from the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/), RNA expression for Glioblastoma Multiforme using Affymetrix HT-HG-U133A (May 14, 2015). Clinical data such as gender, age, histological type, survival and outcome were also downloaded from TCGA data portal. Immune scores and stromal scores were calculated by applying the ESTIMATE algorithm to the downloaded database. For validation, Level 3 gene expression profiles for GBM patients were obtained from the Chinese Glioma Genome Atlas (CGGA) data portal (http://www.cgga.org.cn/), and the RNA sequencing of Diffuse Gliomas using Illumina Hiseq 2000 (Jan 1,2013 and Oct 14 2016). Clinical data of survival and outcome were also downloaded from the CGGA data portal.

22Data analysis was performed using package limma. Fold change > 1.5 and adj. p < 0.05 were set as the cut-offs to screen for differentially expressed genes (DEGs).

23Heatmaps and clustering were generated using an open source web tool ClustVis.

24The protein-protein interaction (PPI) network was retrieved from STRING database and reconstructed via Cytoscape software. Only individual networks with 10 or more nodes were included for further analysis. Networks with fewer than 10 nodes were excluded. The connectivity degree of each node of the network was calculated. Molecular COmplex DEtection (MCODE) was then used to find clusters based on topology to locate densely connected regions.

25Kaplan-Meier plots were generated to illustrate the relationship between patients’ overall survival and gene expression levels of DEGs. The relationship was tested by log-rank test.

26Functional enrichment analysis of DEGs was performed by DAVID (The Database for Annotation, Visualization and Integrated Discovery) to identify GO categories by their biological processes (BP), molecular functions (MF), or cellular components (CC). The DAVID database was also used to perform pathway enrichment analysis with reference from KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways. False discovery rate (FDR) < 0.05 was used as the cut-off.