Mining TCGA database for genes of prognostic value in glioblastoma microenvironment

Results

Immune scores and stromal scores are significantly associated with GBM subtypes

We 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 [6-8], 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.

Figure 1. Immune scores and stromal scores are associated with GBM subtypes and their overall survival. (A) Distribution of immune scores of GBM subtypes. Box-plot shows that there is significant association between GBM subtypes and the level of immune scores (n=417, p<0.001). (B) Distribution of stromal scores of GBM subtypes. Box-plot shows that there is significant association between GBM subtypes and the level of stromal scores (n=417, p<0.001). (C) Distribution of immune scores for IDH1 mutant and IDH1 wildtype GBM cases. Box-plot shows that there is no significant association between IDH1 mutation status and immune scores (n=417, p=0.2758). (D) Distribution of stromal scores for IDH1 mutant and IDH1 wildtype GBM cases. Box-plot shows that there is no significant association between GBM subtypes and the level of stromal scores (n=417, p=0.3077). (E) GBM cases were divided into two groups based on their immune scores: the top half of 209 cases with higher immune scores and the bottom half of 208 cases with lower immune scores. As shown in the Kaplan-Meier survival curve, median survival of the low score group is longer than high score group (442 days vs. 394 days), as indicated by the log-rank test, p value is 0.0537. (F) Similarly, GBM cases were divided into two groups based on their stromal scores: the top half of 209 cases and the bottom half of 208 cases. The median survival of the low score group is longer than the high score group (442 days vs. 422 days), however, it is not statistically different as indicated by the log-rank test p= 0.1262.

Based 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).

To 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.

Comparison of gene expression profile with immune scores and stromal scores in GBM

To 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).

Figure 2. Comparison of gene expression profile with immune scores and stromal scores in GBM. Heatmaps were drawn based on the average linkage method and Pearson distance measurement method. Genes with higher expression are shown in red, lower expression are shown in green, genes with same expression level are in black. (A) Heatmap of the DEGs of immune scores of top half (high score) vs. bottom half (low score). p<0.05, fold change >1.5). (B) Heatmap of the DEGs of stromal scores of top half (high score) vs. bottom half (low score). p<0.05, fold change >1.5). (C, D) Venn diagrams showing the number of commonly upregulated (C) or downregulated (D) DEGs in stromal and immune score groups. (E, F, G) Top 10 GO terms. False discovery rate (FDR) of GO analysis was acquired from DAVID functional annotation tool. p <0.05.

To 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).

Correlation of expression of individual DEGs in overall survival

To 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).

Figure 3. Correlation of expression of individual DEGs in overall survival in TCGA. Kaplan-Meier survival curves were generated for selected DEGs extracted from the comparison of groups of high (red line) and low (blue line) gene expression. p<0.05 in Log-rank test. OS, overall survival in days.

Protein-protein interactions among genes of prognostic value

To 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.

Figure 4. Top 3 PPI networks of IL6, TIMP1, and TLR2 modules. The color of a node in the PPI network reflects the log (FC) value of the Z score of gene expression, and the size of node indicates the number of interacting proteins with the designated protein.

Functional enrichment analysis of genes of prognostic value

Consistent 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.

Figure 5. GO term and KEGG pathway analysis for DEGs significantly associated with overall survival. Top pathways with FDR < 0.05, -log FDR >1.301 are shown: (A) biological process, (B) cellular component, (C) molecular function, and (D) KEGG pathway.

Validation in the CGGA database

To 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.

Figure 6. Validation of correlation of DEGs extracted from TCGA database with overall survival in CGGA cohort. Kaplan-Meier survival curves were generated for selected DEGs extracted from the comparison of groups of high (red line) and low (blue line) gene expression. p<0.05 in Log-rank test. OS, overall survival in days.

Table 1. Genes significant in GBM overall survival identified in both TCGA and CGGA.

Categories	Gene symbols
Cell surface	CD163, CD33, TREM1, THBD
Chemokines		CCL2, CCL5,CCL18, CCL20, CCR5, CXCR4
Complements		C1QA, C1QB, CFH, C3AR1, C5AR1, VSIG4
Tumor necrosis factors		TNFAIP2, TNFAIP6, TNFRSF1B
Interleukins		IL1B, IL6, IL8, IL1R2, IL7R, IL13RA2
Toll-like receptors		TLR1, TLR2
Serpins		SERPING1, SERPINE1
Extracellular matrices	COL1A2, COL5A1, COL6A2,COL6A3, LAMB1, FNDC3B, SRGN
ECM enzymes	TIMP1, CTSS, ADAMTS1, PCOLCE
Cell adhesion molecules	CD44, ICAM1, ITGB2, THBS1
*Genes in bold have not been previously reported for their prognostic value in GBM patients.

Discussion

In 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).

Figure 7. Work flow of the current work.

First, 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 [23-27].

Next, 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 [28-34], indicating poor prognosis.

Finally, 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.

We 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 [35]. As an immunoglobulin supergene family member, ICAM1 is shown to be critically involved in adhesion of cancer cells to ECM [36], in particular, overexpression of ICAM1 correlates with increased tumor malignancy and poor outcome in lung cancer [37], clear cell renal cell carcinoma [38], and mouse GBM model [39].

Significant 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 [6-8,40,41].

The 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 [6]. 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.

In 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.

Acknowledgments

The authors thank Dr. Joanna S. O'Leary for editing the manuscript.

Conflicts of Interest

The authors declare that they have no competing interests.

Funding

This work was supported by the start-up fund from the Vivian L. Smith Department of Neurosurgery, McGovern Medical School, the University of Texas Health Science Center at Houston, Memorial Hermann Foundation-Staman Ogilvie Fund [to Y.L.], the Bentsen Stroke Center [to Y.L.], the TIRR Foundation through Mission Connect [014-115, 017-113 to Y.L.], and the Craig H. Neilsen Foundation [338617 to Y.L.].

References

• 1. Stupp R, Mason WP, van den Bent MJ, Weller M, Fisher B, Taphoorn MJ, Belanger K, Brandes AA, Marosi C, Bogdahn U, Curschmann J, Janzer RC, Ludwin SK, et al, and European Organisation for Research and Treatment of Cancer Brain Tumor and Radiotherapy Groups, and National Cancer Institute of Canada Clinical Trials Group. Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N Engl J Med. 2005; 352:987–96. https://doi.org/10.1056/NEJMoa043330 [PubMed]
• 2. Perry JR, Laperriere N, O’Callaghan CJ, Brandes AA, Menten J, Phillips C, Fay M, Nishikawa R, Cairncross JG, Roa W, Osoba D, Rossiter JP, Sahgal A, et al, and Trial Investigators. Short-Course Radiation plus Temozolomide in Elderly Patients with Glioblastoma. N Engl J Med. 2017; 376:1027–37. https://doi.org/10.1056/NEJMoa1611977 [PubMed]
• 3. Aldape K, Zadeh G, Mansouri S, Reifenberger G, von Deimling A. Glioblastoma: pathology, molecular mechanisms and markers. Acta Neuropathol. 2015; 129:829–48. https://doi.org/10.1007/s00401-015-1432-1 [PubMed]
• 4. Cancer Genome Atlas Research Network. Comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature. 2008; 455:1061–68. https://doi.org/10.1038/nature07385 [PubMed]
• 5. Verhaak RG, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, Miller CR, Ding L, Golub T, Mesirov JP, Alexe G, Lawrence M, O’Kelly M, et al, and Cancer Genome Atlas Research Network. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell. 2010; 17:98–110. https://doi.org/10.1016/j.ccr.2009.12.020 [PubMed]
• 6. Wang Q, Hu B, Hu X, Kim H, Squatrito M, Scarpace L, deCarvalho AC, Lyu S, Li P, Li Y, Barthel F, Cho HJ, Lin YH, et al. Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates with Immunological Changes in the Microenvironment. Cancer Cell. 2017; 32:42–56.e6. https://doi.org/10.1016/j.ccell.2017.06.003 [PubMed]
• 7. Brennan CW, Verhaak RG, McKenna A, Campos B, Noushmehr H, Salama SR, Zheng S, Chakravarty D, Sanborn JZ, Berman SH, Beroukhim R, Bernard B, Wu CJ, et al, and TCGA Research Network. The somatic genomic landscape of glioblastoma. Cell. 2013; 155:462–77. https://doi.org/10.1016/j.cell.2013.09.034 [PubMed]
• 8. Li J, Lu Y, Akbani R, Ju Z, Roebuck PL, Liu W, Yang JY, Broom BM, Verhaak RG, Kane DW, Wakefield C, Weinstein JN, Mills GB, Liang H. TCPA: a resource for cancer functional proteomics data. Nat Methods. 2013; 10:1046–47. https://doi.org/10.1038/nmeth.2650 [PubMed]
• 9. Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, Ohgaki H, Wiestler OD, Kleihues P, Ellison DW. The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary. Acta Neuropathol. 2016; 131:803–20. https://doi.org/10.1007/s00401-016-1545-1 [PubMed]
• 10. Suvà ML, Rheinbay E, Gillespie SM, Patel AP, Wakimoto H, Rabkin SD, Riggi N, Chi AS, Cahill DP, Nahed BV, Curry WT, Martuza RL, Rivera MN, et al. Reconstructing and reprogramming the tumor-propagating potential of glioblastoma stem-like cells. Cell. 2014; 157:580–94. https://doi.org/10.1016/j.cell.2014.02.030 [PubMed]
• 11. Curry JM, Sprandio J, Cognetti D, Luginbuhl A, Bar-ad V, Pribitkin E, Tuluc M. Tumor microenvironment in head and neck squamous cell carcinoma. Semin Oncol. 2014; 41:217–34. https://doi.org/10.1053/j.seminoncol.2014.03.003 [PubMed]
• 12. Cooper LA, Gutman DA, Chisolm C, Appin C, Kong J, Rong Y, Kurc T, Van Meir EG, Saltz JH, Moreno CS, Brat DJ. The tumor microenvironment strongly impacts master transcriptional regulators and gene expression class of glioblastoma. Am J Pathol. 2012; 180:2108–19. https://doi.org/10.1016/j.ajpath.2012.01.040 [PubMed]
• 13. Galon J, Pagès F, Marincola FM, Thurin M, Trinchieri G, Fox BA, Gajewski TF, Ascierto PA. The immune score as a new possible approach for the classification of cancer. J Transl Med. 2012; 10:1. https://doi.org/10.1186/1479-5876-10-1 [PubMed]
• 14. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
• 15. Şenbabaoğlu Y, Gejman RS, Winer AG, Liu M, Van Allen EM, de Velasco G, Miao D, Ostrovnaya I, Drill E, Luna A, Weinhold N, Lee W, Manley BJ, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 2016; 17:231. https://doi.org/10.1186/s13059-016-1092-z [PubMed]
• 16. Winslow S, Lindquist KE, Edsjö A, Larsson C. The expression pattern of matrix-producing tumor stroma is of prognostic importance in breast cancer. BMC Cancer. 2016; 16:841. https://doi.org/10.1186/s12885-016-2864-2 [PubMed]
• 17. Hanahan D, Weinberg RA. The hallmarks of cancer. Cell. 2000; 100:57–70. https://doi.org/10.1016/S0092-8674(00)81683-9 [PubMed]
• 18. Hanahan D, Coussens LM. Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell. 2012; 21:309–22. https://doi.org/10.1016/j.ccr.2012.02.022 [PubMed]
• 19. Carter SL, Cibulskis K, Helman E, McKenna A, Shen H, Zack T, Laird PW, Onofrio RC, Winckler W, Weir BA, Beroukhim R, Pellman D, Levine DA, et al. Absolute quantification of somatic DNA alterations in human cancer. Nat Biotechnol. 2012; 30:413–21. https://doi.org/10.1038/nbt.2203 [PubMed]
• 20. Shah N, Wang P, Wongvipat J, Karthaus WR, Abida W, Armenia J, Rockowitz S, Drier Y, Bernstein BE, Long HW, Freedman ML, Arora VK, Zheng D, Sawyers CL. Regulation of the glucocorticoid receptor via a BET-dependent enhancer drives antiandrogen resistance in prostate cancer. eLife. 2017; 6:e27861. https://doi.org/10.7554/eLife.27861 [PubMed]
• 21. Priedigkeit N, Watters RJ, Lucas PC, Basudan A, Bhargava R, Horne W, Kolls JK, Fang Z, Rosenzweig MQ, Brufsky AM, Weiss KR, Oesterreich S, Lee AV. Exome-capture RNA sequencing of decade-old breast cancers and matched decalcified bone metastases. JCI Insight. 2017; 2:95703. https://doi.org/10.1172/jci.insight.95703 [PubMed]
• 22. Alonso MH, Aussó S, Lopez-Doriga A, Cordero D, Guinó E, Solé X, Barenys M, de Oca J, Capella G, Salazar R, Sanz-Pamplona R, Moreno V. Comprehensive analysis of copy number aberrations in microsatellite stable colon cancer in view of stromal component. Br J Cancer. 2017; 117:421–31. https://doi.org/10.1038/bjc.2017.208 [PubMed]
• 23. Silver DJ, Sinyuk M, Vogelbaum MA, Ahluwalia MS, Lathia JD. The intersection of cancer, cancer stem cells, and the immune system: therapeutic opportunities. Neuro-oncol. 2016; 18:153–59. https://doi.org/10.1093/neuonc/nov157 [PubMed]
• 24. Lu P, Takai K, Weaver VM, Werb Z. Extracellular matrix degradation and remodeling in development and disease. Cold Spring Harb Perspect Biol. 2011; 3:a005058. https://doi.org/10.1101/cshperspect.a005058 [PubMed]
• 25. Bonnans C, Chou J, Werb Z. Remodelling the extracellular matrix in development and disease. Nat Rev Mol Cell Biol. 2014; 15:786–801. https://doi.org/10.1038/nrm3904 [PubMed]
• 26. Pickup MW, Mouw JK, Weaver VM. The extracellular matrix modulates the hallmarks of cancer. EMBO Rep. 2014; 15:1243–53. https://doi.org/10.15252/embr.201439246 [PubMed]
• 27. Peranzoni E, Rivas-Caicedo A, Bougherara H, Salmon H, Donnadieu E. Positive and negative influence of the matrix architecture on antitumor immune surveillance. Cell Mol Life Sci. 2013; 70:4431–48. https://doi.org/10.1007/s00018-013-1339-8 [PubMed]
• 28. Yeung YT, McDonald KL, Grewal T, Munoz L. Interleukins in glioblastoma pathophysiology: implications for therapy. Br J Pharmacol. 2013; 168:591–606. https://doi.org/10.1111/bph.12008 [PubMed]
• 29. Tchirkov A, Khalil T, Chautard E, Mokhtari K, Véronèse L, Irthum B, Vago P, Kémény JL, Verrelle P. Interleukin-6 gene amplification and shortened survival in glioblastoma patients. Br J Cancer. 2007; 96:474–76. https://doi.org/10.1038/sj.bjc.6603586 [PubMed]
• 30. Wang XB, Tian XY, Li Y, Li B, Li Z. Elevated expression of macrophage migration inhibitory factor correlates with tumor recurrence and poor prognosis of patients with gliomas. J Neurooncol. 2012; 106:43–51. https://doi.org/10.1007/s11060-011-0640-3 [PubMed]
• 31. Cuny E, Loiseau H, Penchet G, Ellie E, Arsaut J, Vital A, Vincendeau P, Demotes-Mainard J. Association of elevated glial expression of interleukin-1beta with improved survival in patients with glioblastomas multiforme. J Neurosurg. 2002; 96:294–301. https://doi.org/10.3171/jns.2002.96.2.0294 [PubMed]
• 32. Vinnakota K, Hu F, Ku MC, Georgieva PB, Szulzewsky F, Pohlmann A, Waiczies S, Waiczies H, Niendorf T, Lehnardt S, Hanisch UK, Synowitz M, Markovic D, et al. Toll-like receptor 2 mediates microglia/brain macrophage MT1-MMP expression and glioma expansion. Neuro-oncol. 2013; 15:1457–68. https://doi.org/10.1093/neuonc/not115 [PubMed]
• 33. Cai J, Zhang W, Yang P, Wang Y, Li M, Zhang C, Wang Z, Hu H, Liu Y, Li Q, Wen J, Sun B, Wang X, et al. Identification of a 6-cytokine prognostic signature in patients with primary glioblastoma harboring M2 microglia/macrophage phenotype relevance. PLoS One. 2015; 10:e0126022. https://doi.org/10.1371/journal.pone.0126022 [PubMed]
• 34. Pan Y, Smithson LJ, Ma Y, Hambardzumyan D, Gutmann DH. Ccl5 establishes an autocrine high-grade glioma growth regulatory circuit critical for mesenchymal glioblastoma survival. Oncotarget. 2017; 8:32977–89. https://doi.org/10.18632/oncotarget.16516 [PubMed]
• 35. Rajaraman P, Brenner AV, Butler MA, Wang SS, Pfeiffer RM, Ruder AM, Linet MS, Yeager M, Wang Z, Orr N, Fine HA, Kwon D, Thomas G, et al. Common variation in genes related to innate immunity and risk of adult glioma. Cancer Epidemiol Biomarkers Prev. 2009; 18:1651–58. https://doi.org/10.1158/1055-9965.EPI-08-1041 [PubMed]
• 36. Nakada M, Nakada S, Demuth T, Tran NL, Hoelzinger DB, Berens ME. Molecular targets of glioma invasion. Cell Mol Life Sci. 2007; 64:458–78. https://doi.org/10.1007/s00018-007-6342-5 [PubMed]
• 37. Kotteas EA, Boulas P, Gkiozos I, Tsagkouli S, Tsoukalas G, Syrigos KN. The intercellular cell adhesion molecule-1 (icam-1) in lung cancer: implications for disease progression and prognosis. Anticancer Res. 2014; 34:4665–72. [PubMed]
• 38. Roland CL, Harken AH, Sarr MG, Barnett CCJr. ICAM-1 expression determines malignant potential of cancer. Surgery. 2007; 141:705–07. https://doi.org/10.1016/j.surg.2007.01.016 [PubMed]
• 39. Piao Y, Henry V, Tiao N, Park SY, Martinez-Ledesma J, Dong JW, Balasubramaniyan V, de Groot JF. Targeting intercellular adhesion molecule-1 prolongs survival in mice bearing bevacizumab-resistant glioblastoma. Oncotarget. 2017; 8:96970–83. https://doi.org/10.18632/oncotarget.18859 [PubMed]
• 40. Brat DJ, Verhaak RG, Aldape KD, Yung WK, Salama SR, Cooper LA, Rheinbay E, Miller CR, Vitucci M, Morozova O, Robertson AG, Noushmehr H, Laird PW, et al, and Cancer Genome Atlas Research Network. Comprehensive, Integrative Genomic Analysis of Diffuse Lower-Grade Gliomas. N Engl J Med. 2015; 372:2481–98. https://doi.org/10.1056/NEJMoa1402121 [PubMed]
• 41. Ceccarelli M, Barthel FP, Malta TM, Sabedot TS, Salama SR, Murray BA, Morozova O, Newton Y, Radenbaugh A, Pagnotta SM, Anjum S, Wang J, Manyam G, et al, and TCGA Research Network. Molecular Profiling Reveals Biologically Discrete Subsets and Pathways of Progression in Diffuse Glioma. Cell. 2016; 164:550–63. https://doi.org/10.1016/j.cell.2015.12.028 [PubMed]
• 42. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
• 43. Metsalu T, Vilo J. ClustVis: a web tool for visualizing clustering of multivariate data using Principal Component Analysis and heatmap. Nucleic Acids Res. 2015; 43:W566-70. https://doi.org/10.1093/nar/gkv468 [PubMed]
• 44. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, Kuhn M, Bork P, Jensen LJ, von Mering C. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015; 43:D447–52. https://doi.org/10.1093/nar/gku1003 [PubMed]
• 45. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
• 46. Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009; 4:44–57. https://doi.org/10.1038/nprot.2008.211 [PubMed]