Bioinformatic Profiling of Prognosis-Related Genes in Malignant Glioma Microenvironment.

BACKGROUND Gliomas are the most common primary tumors of the brain and spinal cord. The tumor microenvironment (TME) is the cellular environment in which tumors exist. This study aimed to identify the role of the TME and the effects of genes involved in the TME of malignant glioma. MATERIAL AND METHODS The ESTIMATE algorithms in the R package were used to calculate the immune and stromal scores of samples in the TCGA and GSE4290 datasets. The associations of stromal and immune scores with clinicopathological characteristics and overall survival of malignant glioma patients were assessed by analysis of variance and Kaplan-Meier analysis. Differentially expressed genes (DEGs) were obtained through the median immune and stromal score using the R package "limma". Functional enrichment analysis and the PPI network MCODE were used to analyze DEGs. RESULTS Increased immune and stromal scores were closely related with advanced glioma grade and poor prognosis (all P<0.01). In total, 558 DEGs were found and most were related to tumor prognosis. Functional enrichment analysis showed that DEGs were associated with cell-matrix regulation and immune response. Four hub modules related to tumor angiogenesis, collagen formation, and immune response were found and analyzed. Previously overlooked microenvironment-related genes such as LAMB1, FN1, ACTN1, TRIM, SERPINH1, CYBA, LAIR1, and LILRB2 showed prognostic values in malignant glioma patients. CONCLUSIONS The glioma stromal/immune scores are closely related to glioma grade, histology, and survival time. Some glioma microenvironment-related genes including LAMB1, FN1, ACTN1, TRIM6, SERPINH1, CYBA, LAIR1, and LILRB2 show prognostic values in malignant gliomas and serve as potential biomarkers.

Background

Gliomas which originate from astrocyte or oligodendrocyte precursor cells are the most common primary tumors in the brain and spinal cord [1]. Gliomas comprise about 30% of all brain tumors and 80% of malignant brain tumors [2]. Gliomas are divided into 4 grades (World Health Organization [WHO] grades I to IV) according to the biological behavior and malignancy of the tumor [3]. WHO grade II–IV gliomas, which are referred to as malignant gliomas, show an invasive growth mode and are difficult to cure by surgical treatment, radiotherapy, chemotherapy, or other means [4]. Although great strides have been made in the current comprehensive treatment of malignant gliomas, the overall prognosis of patients with malignant glioma, especially grade IV glioblastoma multifore (GBM), has not significantly improved. The overall survival of patients with GBM is only about 14.2 months [5]. At present, surgery cannot achieve total tumor resection, and the combination of surgery with postoperative radiotherapy and chemotherapy is still the preferred option [6]. Therefore, there is an urgent need to explore and find new treatment methods for malignant gliomas to improve the prognosis and quality of life of patients [7].

In recent years, the relationship between a tumor and its surrounding tissue has been widely discussed. Tumor surrounding tissue which is called the “tumor microenvironment” (TME) is made up of extracellular components including cytokines, growth factors, hormones, and extracellular matrix, and different types of cells including endothelial cells, fibroblasts, and immune cells [8].The extracellular matrix exists in all tissues and constitutes the non-cellular component of the microenvironment and functions as a physical scaffold and a source of biochemical signals. It maintains a close relationship with intracellular biochemical and biomechanical processes and strongly influences the biochemistry and biomechanics processes of the tissue. Cooperation between cellular and extracellular matrix factors regulates cellular fate, tissue morphology, and organogenesis [9]. The extracellular matrix constitutes about 20% of brain volume, varying in composition among different brain regions [10]. Like that of most solid tumors, the glioma TME contains other cell types in addition to neoplastic tumor cells, in particular vascular cells and immune cells including CD4+ and CD8+ T lymphocytes, T regulatory lymphocytes, glioma associated microglia/macrophages, dendritic cells, myeloid-derived suppressor cells, and natural killer cells [11]. These cells, which are attracted by soluble factors secreted by neurocytes, normally regulate tissue homeostasis, immune surveillance, and wound healing. But immune cells in glioma patients are often disabled, failing to activate effector T cells and durable anti-tumor immune responses [12]. In view of the defect of immune response in malignant gliomas, novel therapies for immune responses are at present being developed, and the glioma TME has received increasing attention as a potential therapeutic target [13-15].

Previously, pathologists estimated the composition of tumors using visual assessments, while with the development of bioinformation, a variety of algorithms have been developed to investigate the interactions between tumor cells and the TME in tumor tissues. Yoshihara et al. developed an ESTIMATE algorithm that can calculate the molecular markers of immune cells and stromal cells in tumor tissues to get the immune and stromal ESTIMATE score of the tumor to predict the TME [16]. A high immune score indicates there are more immune cell components, while a high stromal score indicates there are more stromal components. The ESTIMATE score is the sum of the immune and stromal score. The higher the ESTIMATE score, the more immune and stromal components the tumor has, and the lower the purity of tumor. Based on ESTIMATE algorithms, researchers have explored prognostic assessment and genetic changes in various tumors [17-19], and it has been experimentally confirmed that the algorithms are consistent with the actual situation of the tissue surrounding the tumor. This allows us to calculate the immune and stromal scores to predict and assess the TME in gliomas.

In this study, we performed datamining in the Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), and Chinese Glioma Genome Atlas (CGGA) to download and investigate the expression profile and clinical information of malignant gliomas (low grade gliomas and GBM) to determine the prognosis-related genes in malignant glioma TME.

Material and Methods

Data processing

The gene expression profiles of GSE4290 based on the GPL570 platform ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array) were downloaded from the GEO database () [20]. The dataset included 180 samples, 23 controls, and 157 grade II–IV gliomas. A total of 153 glioma samples with histological data were selected to calculate immune and stromal ESTIMATE scores using the R package “estimate”. The TCGA Level 3 data of low grade glioma and GBM were downloaded from the TCGA Data Coordinating Center () [21]. The R package “limma” was used to normalize the raw data, and the package “estimate” was used to calculate the immune, stromal, and ESTIMATE scores. For validation, 325 malignant glioma samples including RNA sequencing and clinical data was downloaded from the CGGA part C data () [22,23].

Identification of differentially expressed genes

The differentially expressed genes (DEGs) of GSE4290 and TCGA malignant glioma data were calculated using the R package “limma” by the median immune and stromal score separately. The cut-off criterion was set as |log2FC| >1.0 and FDR <0.05 to determine the DEGs in the limma package. We used the intersection of the highly expressed immune and stromal genes in GSE4290 and TCGA to get the up-regulated genes in malignant glioma. Then same was done in the low expression group.

Kaplan-Meier analysis and functional enrichment of DEGs

Kaplan-Meier (K-M) curves were generated to illustrate the relationship between patients’ overall survival and gene expression levels of the DEGs by the R package “survival” in the TCGA and CGGA datasets. We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the DEGs with the R package “enrichplot”. We used FDR <0.05 as the cut-off criteria to screen out the enriched GO terms and KEGG pathways.

Protein–protein network construction

We then used the Search Tool for the Retrieval of Interacting Genes (STRING; ) (version 10.0) online database to predict the protein-protein interaction (PPI) network of DEGs and analyze the degree of interactions between proteins with a combined score >0.9 [24]. The network was then visualized by using Cytoscape software (version 3.7.1).

Identification of hub modules

The hub modules were identified by Molecular Complex Detection (MCODE) (version1.5.1), a plug-in of Cytoscape. MCODE could identify the most significant module in the PPI networks by clustering the given network based on topology to find densely connected regions [25]. We used degree cutoff=2, k-score=2, node score cutoff=0.2, and max depth=100 as the criteria to seek the hub modules. Then GO analysis was performed in every module.

Statistical analyses

All statistical analyses were performed in R version 3.3.6 (http//). The associations of stromal and immune scores with clinicopathological characteristics and overall survival of malignant glioma patients were assessed by analysis of variance and K-M analysis. Results were considered statistically significant for a 2-tailed P value <0.05.

Results

Increased immune and stromal scores closely correlated with higher glioma grade and poor prognosis

In the study, we obtained the expression profile of 665 malignant gliomas in TCGA which had complete clinical data. We first converted the probe ID to gene name, and then averaged the repeated gene expression. Then, we processed the ESTIMATE algorithm and got the stromal score, immune score, and ESTIMATE score (Supplementary Table 1). K-M analysis was done to reveal the correlation between the scores and survival by dividing 665 patients into high and low score groups with the median. Survival curves indicated that high immune score significantly correlated with poor overall survival (P<0.001; Figure 1A). The same trend occurred in the high stromal (P<0.001) and ESTIMATE (P<0.001) score groups (Figure 1B, 1C).

Figure 1

Association between immune, stromal, and Estimate score and survival in TCGA. (A) Survival curves indicate that elevated immune score significantly correlated with poor overall survival (P<0.001). (B) Elevated stromal score significantly associated with poor overall survival (P<0.001). (C) Elevated ESTIMATE score significantly associated with poor overall survival (P<0.001).

Then, we evaluated the correlation between the scores and clinicopathological indicators. Immune and stromal scores gradually increased with the increase of glioma grade (both P<0.001), and grade IV gliomas had the highest score and poorest overall survival (Figure 2A, 2B). Based on the histology, there were significant differences among different histologic gliomas (both P<0.001), among which GBM had the highest score (Figure 2C, 2D). Patients with IDH/1-mutated gliomas had a lower stromal score, while immune score had no significant relationship with IDH/1 status (Figure 2I, 2J).

Figure 2

Association between immune and stromal scores and clinicopathological characteristics. (A, C) With the increase of glioma grade, the immune and stroma score increased in TCGA. (B, D) Glioblastoma has the highest immune and stroma score in TCGA. (E, G) With the increase of glioma grade, the immune and stroma scores increased in GSE4290. (F, H) Glioblastoma has the highest immune and stroma scores in GSE4290. (I, J) IDH/1mutant patients have higher stromal scores but not higher immune scores.

Using the same algorithm, we scored 153 grade II–IV gliomas in GSE4290 (Supplementary Table 2). Due to the lack of survival data, we used it only to verify the relationship between scores and clinicopathological indicators. The results in GSE4290 were very similar to those in TCGA whereby immune and stromal scores increased with the advanced glioma grade, and GBM had the highest score of all histological types (all P<0.001; Figure 2E–2H).

Identification of DEGs with immune and stromal scores in malignant glioma

To explore DEG profiles, we compared 665 cases obtained in the TCGA and 153 cases in the GSE4290 dataset by immune and stromal score. In the TCGA dataset, based on the cut-off criteria (P<0.05 and |log2FC| >1.0), a total of 4033 DEGs, including 2829 up-regulated and 1204 down-regulated genes, were identified between high and low immune scores; 4629 DEGs, including 3392 up-regulated and 1237 down-regulated genes, were identified between high and low stromal scores. In GSE4290, 491 up-regulated and 375 down-regulated genes were identified by immune scores, and 593 up-regulated and 357 down-regulated genes were identified by stromal scores.

Using the Venn diagram to take the intersection of the significantly high expressed immune and stromal genes in GSE4290 and TCGA, we obtained 362 up-regulated and 196 down-regulated genes both in immune and stromal score in the TCGA and GSE4290 datasets (Figure 3 and Supplementary Table 3). These 558 genes were considered DEGs and extracted for further analysis.

Figure 3

Venn algorithm shows differentially expressed genes by immune and stromal scores. (A) A total of 362 DEGs were commonly up-regulated by the immune and stroma scores in TCGA and GSE4290. (B) A total of 196 genes were commonly down-regulated by the immune and stroma scores in TCGA and GSE4290.

Functional enrichment analysis of DEGs

To explore the function of DEGs, we performed GO enrichment analysis. The DEGs were categorized into 3 functional groups in GO terms: biological process, cellular component, and molecular function (Figure 4A and Supplementary Table 4). DEGs in the biological process group were enriched mainly in neutrophil activation (GO: 0042119, P=4.19E-26), neutrophil degranulation (GO: 0043312, P=1.06E-23), neutrophil activation involved in immune response (GO: 0002283, P=1.48E-23), neutrophil mediated immunity (GO: 0002446, P=4.88601825321153E-23), and extracellular matrix organization (GO: 0030198, P=2.214E-18). The genes in the cellular component group were significantly enriched in extracellular matrix (GO: 0031012, P=8.18E-32), collagen-containing extracellular matrix (GO: 0062023, P=8.92E-30), cytoplasmic vesicle lumen (GO: 0060205, P=3.54E-16), vesicle lumen (GO: 0031983, P=3.94E-16), and collagen trimer (GO: 0005581, P=1.55E-15). The molecular function group were enriched in extracellular matrix structural constituent (GO: 0005201, P=1.23E-16), IgG binding (GO: 0019864, P=7.55E-11), extracellular matrix structural constituent conferring tensile strength (GO: GO: 0030020, P=5.00E-09), toll-like receptor binding (GO: 0035325, P=5.12E-09), and immunoglobulin binding (GO: 0019865, P=8.05E-09). KEGG pathway analysis showed that the top enriched terms were phagosome (hsa04145, P=2.54E-18), staphylococcus aureus infection (hsa05150, P=4.11E-16), leishmaniasis (hsa05140, P=2.81E-14), complement and coagulation cascades (hsa04610, 6.57E-14), and tuberculosis (hsa05152, P=6.09E-12) based on P value (Figure 4B and Supplementary Table 5). The results suggest that these genes were closely related to the cell matrix regulation and immune response.

Figure 4

Functional enrichment of DEGs. (A) Gene ontology analysis. (B) KEGG pathway analysis.

Survival analysis of DEGs in TCGA and validation dataset

K-M analysis was used to assess the potential roles of DEGs in overall survival in the TCGA dataset malignant glioma patients (Supplementary Table 6). A total of 543 DEGs were shown to significantly predict overall survival in log-rank tests among the 558 DEGs (P<0.05, selected genes are shown in Figure 5). More convincingly, when we verified the result in the CGGA dataset, 541 of the 558 genes were found in CGGA, 539 of which were confirmed to be significantly associated with prognosis prediction (P<0.05, Supplementary Table 7, and selected genes are shown in Figure 6). The above results show that most of the DEGs we obtained are related to tumor prognosis, which warrants further study.

Figure 5

Overall survival of 9 selected DEGs in TCGA based on Kaplan-Meier plotter. The patients were stratified into a high-level group and a low-level group according to median expression. (A) ACTN1. (B) CYBA. (C) FN1. (D) LAIR1. (E) LAMB1. (F) LILRB2. (G) SERPINH1. (H) SP100. (I) TRIM6.

Figure 6

Overall survival of 9 selected DEGs in CGGA based on Kaplan-Meier plotter. The patients were stratified into a high-level group and a low-level group according to median expression. (A) ACTN1. (B) CYBA. (C) FN1. (D) LAIR1. (E) LAMB1. (F) LILRB2. (G) SERPINH1. (H) SP100. (I) TRIM6.

PPI network construction and identification of hub modules

To better understand the interactions among the DEGs, a PPI network was obtained from the STRING database (Supplementary Figure 1). This network was made up of 278 nodes and 1017 edges. A total of 16 modules were screened out by MCODE. The 4 most significant modules were selected for further analysis (Figure 7 and Supplementary Table 8). The GO analysis results of each module are shown in Figure 8.

Figure 7

(A–D) Top 4 PPI network modules. Red represents the up-regulation of gene expression, while blue represents the down-regulation. Node size indicates the number of interacting proteins with the other DEGs in the original PPI network, and the degree of thickness reflects the combined_score in the STRING dataset.

Figure 8

GO analysis of genes in the 4 significant modules. (A) Module 1. (B) Module 2. (C) Module 3. (D) Module 4.

Discussion

TME was first proposed by Lord in 1979, and a relatively complete theoretical system has since been established [26]. Current researches focus primarily on microenvironment-mediated therapeutic resistance. Soluble factor, cell adhesion, and immune response-mediated drug resistance were observed in tumor treatment failure [27]. Therefore, therapies targeting stroma-derived paracrine factors, tumor cell–extracellular matrix interactions, and immune cell therapy have been developed and clinically tested.

In recent years, due to slow progress in the traditional treatment of gliomas, researchers turned their attention to molecular therapy and immunotherapy. However, due to the obstruction of the blood–brain barrier and the resistance to the TME of the central system, immunotherapy has not made significant progress [28]. In recognition of the important role of TME in treatment resistance, it is believed that blocking TME pathways would prevent or delay the acquisition of tumor progress. Therefore, increasing attention has been given to the glioma microenvironment as a potential therapeutic target.

The development of bioinformatics has greatly changed biomedicine research through the development of microarray sequencing and the establishment of the open database. Researchers increasingly explore new targets for malignant gliomas and use statistical algorithms for external validation [29].

Currently, by extracting the expression profile from TCGA and GSE4290, we calculated stromal, immune, and ESTIMATE scores using the ESTIMATE algorithm. K-M analysis indicated that high immune and stromal scores significantly correlated with poor overall survival. Immune and stromal scores gradually increased with an increase in WHO grade, and grade IV GBM had the highest score and poorest overall survival. Then, we obtained 558 DEGs by high and low immune and stromal scores in TCGA and GSE4290 and most of the DEGs were related to tumor prognosis in TCGA and CGGA. GO and KEGG analysis showed that DEGs were closely related to cell matrix regulation and immune response. A total of 16 significant modules were identified using MCODE, which was based on topology to find densely connected regions in Cytoscape. These modules are composed of genes with significant functions, and each module represents 1 or more physiological processes correlated with the glioma TME. From these modules, we explored the mode of TME to identify the most critical gene components. The 4 most significant modules are highlighted below.

In module 1, it can be seen that C3 and VEGFA are in the center of the module. VEGFA is a member of vascular endothelial growth factor (VEGF), which is often considered to be the prototypical angiogenic molecule in malignancy. Although tumor cells are major sources of VEGF, stromal cells within the TME are additional contributors [30]. Most of the left side of the module is composed of C-C motif chemokine-related genes. Through GO analysis, we found that these genes are located primarily in the external side of the plasma membrane, participating in chemokine activity and G-protein-combined receptor binding. Their functions are mainly leukocyte chemotaxis and regulation. In human colon cancer, CXCR4 is overexpressed in the chemoresistant tumor cells. Stromal cells from lymph nodes promote resistance to 5-fluorouracil and oxaliplatin through a SDF1/CXCR4 dependent mechanism [31]. Many other genes are reported to be related to gliomas, for example FPR1, SSTR1, and ANXA1 [32-34]. The gene right sides are composed of vesicle structures, such as platelet alpha granule, endoplasmic reticulum, and secret granule, which participate in the functions of platelet expansion and post translational protein modification. We are particularly interested in LAMB1, FN1, and ACTN1, which have been reported to play important roles in other tumors, but whose roles are still unclear in gliomas [35-37]. LAMB1, one of extracellular matrix glycoproteins, is the major non-collagenous constituent of basement membranes and has been implicated in a wide variety of biological processes including cell adhesion, differentiation, migration, signaling, neurite outgrowth, and metastasis. Survival analysis in TCGA and CGGA showed that high expression of LAMB1 had significantly lower survival time in malignant gliomas, which may be related to its role in promoting tumor cell migration. FN1 encodes fibronectin, a glycoprotein present in a soluble dimeric form in plasma, and in a dimeric or multimeric form at the cell surface and extracellular matrix. Relevant studies have shown that its function may be to gather around neovascularization, surround cancer cells, support tumor growth, and promote invasion [38,39]. ACTN1 belong to the spectrin gene superfamily, which represents a diverse group of cytoskeletal proteins [40]. In general, module 1 is mainly related to tumor angiogenesis.

Module 2 includes primarily the 2 Fc fragment of IgG receptor genes, 3 guanylate-binding proteins, and some HLA-related genes. They are mainly involved in the formation of clathrin-coated endocytic vesicles and the MHC protein complex in cells. They have the function of binding amine, peptide, and antigen. They participate in response to interference gamma and antigen processing and presentation in the body, which are important parts of the immune response process. Among the remaining genes, CD44 is thought to facilitate glioma interaction with the TME to promote malignancy; clinical trials of CD44 targeted therapy in CD44 positive solid tumors are under way [41]. GILT is a lysosomal thiol reductase. Elevated GILT is positively associated with glioma progression, which is consistent with our analysis [42]. However, it is strange that previous studies reported that SP100 expression can reduce malignancy of brain tumors. But in this analysis, we found that SP100 high expression group has a significantly low survival rate in the TCGA and CGGA cohort [43]. Further study is needed to explore the specific role of SP100 in glioma. TRIM6 belongs to a large family of E3-Ub ligases, which are characterized by the presence of RING, B-box, and coiled-coil domains and have been implicated in innate immunity [44]. No studies to date have exposed its role in gliomas; however our study shows that this molecule may also play an important role in the glioma TME.

Module 3 is primarily composed of collagen genes, which make up the extracellular matrix structure and participate in tissue regulation and collagen fabric organization. In solid tumors, collagen provides mechanical strength and promotes the binding properties of cytokines. A well-organized and interconnected collagen network reduces the amount of drug penetration. In addition, drugs are isolated by binding with extracellular matrix components, thus inhibiting drug penetration into deeper tumor areas [45]. SERPINH1, also named HSP47, is a collagen-specific molecular chaperone that localizes in the endoplasmic reticulum. SERPINH1 is indispensable for molecular maturation of collagen and is reported to enhance glioma tumor growth and invasion [46].

The genes in module 4 form a granular membrane and participate in neutrophil degradation, activation, and immune response. Among them, PLAU, CLEC5A, CD36, and CD93 have been reported to affect the generation and development of gliomas [47-49]. VAMP8 is thought to promote the proliferation of glioma cells and lead to temozolomide resistance in human glioma cells [50]. Other molecules, such as BST1, CYBA, LAIR1, and LILRB2, play an important role in the process of TME, but have not been reported in malignant glioma.

Conclusions

In conclusion, through the ESTIMATE algorithm, differentially expressed analysis, and PPI network analysis, we obtained 4 hub modules related to the glioma TME. By analyzing the modules in detail, we identified previously neglected genes in the glioma microenvironment, namely LAMB1, FN1, ACTN1, TRIM, SERPINH1, CYBA, LAIR1, and LILRB2. Further study on these genes and their roles in gliomas is needed.

Protein-protein network of commonly up-regulated and down-regulated DEGs in both stromal and immune score groups, generated with the STRING database.