An Autophagy-Related Gene Signature Associated With Clinical Prognosis and Immune Microenvironment in Gliomas

Glioma is one of the leading causes of death from cancer, and autophagy-related genes (ARGs) play an important role in glioma occurrence, progression, and treatment. In this study, the gene expression profiles and clinical data of glioma patients were obtained from The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA), respectively. ARGs were obtained from the Human Autophagy Database. We analyzed the expression of the ARGs in glioma and found that 73 ARGs were differentially expressed in tumor and normal tissues. Univariate Cox regression analysis was used to identify prognostic differentially expressed ARGs (PDEARGs). Least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression analyses were performed on the PDEARGs to determine the risk genes; and BRIC5, NFE2L2, GABARAP, IKBKE, BID, MAPK3, FKBP1B, MAPK8IP1, PRKCQ, CX3CL1, NPC1, HSP90AB1, DAPK2, SUPT20H, and PTEN were selected to establish a prognostic risk score model for TCGA and CGGA cohorts. This model accurately stratified patients with different survival outcomes, and the autophagy-related signature was also appraised as being an independent prognostic factor. We also constructed a prognostic nomogram using risk score, age, gender, WHO grade, isocitrate dehydrogenase (IDH) mutation status, and 1p/19q co-deletion status; and the calibration plots showed excellent prognostic performance. Finally, Pearson correlation analysis suggested that the ARG signature also played an essential role in the tumor immune microenvironment. In summary, we constructed and verified a novel autophagy-related signature that was tightly associated with the tumor immune microenvironment and could serve as an independent prognostic biomarker in gliomas.

Introduction

Glioma is one of the most common primary brain tumors in adults, accounting for about 15% of all brain tumors and more than 48% of malignant brain tumors (1). The prognosis of patients with glioma is poor, despite the use of multi-mode therapy, including neurosurgery, radiotherapy, chemotherapy, targeted therapy, and immunotherapy, due to the high recurrence and mortality rates (2, 3). For glioma patients, it still appears difficult to predict and extend survival, and the identification of novel prognostic biomarkers for glioma is therefore imperative. Autophagy is the process of degradation and circulation of proteins and intracellular components during starvation or stress (4). It occurs at the basal level of all cells and is fundamental to cell homeostasis (5, 6). Autophagy often plays a key role in the pathogenesis of diseases, such as type 2 diabetes, cardiovascular disease, neurological disease, and cancer (7). Autophagy also has a dual function in tumors, eliminating damaged proteins and organelles in the early stages of the tumor, reducing cell damage, and thus inhibiting tumor development (8). However, once tumors have formed, autophagy via HIF1α/AMPK can provide nutrients to promote tumor cell growth when nutrients are limited (9). In most cases, autophagy is thought to inhibit early tumorigenesis but promote the development of established tumors.

To illustrate the importance of autophagy, there is increasing evidence of the role of autophagy-related genes (ARGs) in the development of cancer. High expression of Atg10 is associated with poor prognosis for lung cancer (10), and autophagy gene Atg7 knockout leads to the progression of lung cancer to eosinophilic cell tumors, inhibiting tumor proliferation (11). Autophagy has been negatively correlated with mTOR substrate phosphorylation, and recent studies have shown that high expression of LC3/beclin 1 is associated with poor prognosis in glioma patients (12). In addition, autophagy-related proteins, such as cathepsin D and LAMP2, are overexpressed in glioma and may also be potential targets for prognosis and prediction (13). These findings confirm the role of autophagy in cancer and suggest that ARGs may serve as prognostic markers for glioma.

Although some previous studies have explored the relationship between ARGs and prognosis in glioma patients, those studies have focused on the function of a single gene. Very few studies have used gene expression profiles to examine the relationship between multiple ARGs and glioma prognosis. Therefore, in this study, we develop a reliable prognostic model based on ARGs and explore the clinical application of the ARG signatures in glioma patients.

Materials and Methods

Acquisition of Glioma Datasets

The RNA-Seq and somatic mutation data of 703 patients were downloaded from The Cancer Genome Atlas (TCGA) data portal (14)¹ and used as a training cohort. From the Chinese Glioma Genome Atlas (CGGA) dataset², we collected the RNA-Seq data of 1,018 samples as a validation cohort (15, 16). We compared transcriptome data with the patients’ clinical information based on the patient ID numbers. If a patient’s ID number did not match, we removed it. Thus, 665 patients with complete overall survival (OS) information and gene expression profiles were obtained from TCGA database, while 929 patients with complete OS information and gene expression profiles were obtained from the CGGA database. The 232 ARGs were obtained from the human autophagy HADb database³ (17). Immune infiltrate data for glioma patients were derived through a deconvolution algorithm based on a normalized gene expression profile (CIBERSORT), which was used to characterize the infiltrating components of the immune cells in 22 complex samples (18, 19).

Identification of Differentially Expressed Autophagy-Related Genes

The R package “edgeR” was used to screen differentially expressed genes (DEGs) in glioma and normal neural tissue samples (20). A Wilcoxon signed-rank test was used to screen DEGs according to the cutoff values’ false discovery rate (FDR) <0.05 and |log2 FC| > 1. The intersection of DEGs and ARGs was considered to be the set of significant differentially expressed ARGs (DEARGs).

Construction and Verification of a Prognostic Risk Score Model Based on Autophagy-Related Genes

We chose the data in TCGA dataset as the training cohort, which was used to construct the Cox proportional hazard regression model, and the data in the CGGA dataset was used as the validation cohort, which was used to validate the model’s performance. First, univariate Cox analysis was used to identify possible prognostic DEARGs (PDEARGs) by using the “survival” package with a P-value of <0.05 (21). Second, the potential risk genes were selected, and the overfitting genes were eliminated by least absolute shrinkage and selection operator (LASSO) regression using the “glmnet” R package (22). Third, Cox proportional hazard regression was used to construct the prognostic risk model by using the “glmnet” R package (22). The following formula was used to establish the prognostic risk model:

$$riskscore=\sum _{j=1}^n(Coe{f}_j\times X_j),$$

where Coef is the multivariate Cox regression analysis coefficient of the ARGs and X is the relative expression level of each ARG. The patients were divided into high-risk and low-risk groups, with the median risk score used as the cutoff value and high-risk scores suggesting that the prognosis of the patients was poor.

Statistical Analysis

R software (3.6.1) was used for all statistical analyses. Pearson correlation coefficients were used to evaluate the rank correlations among the variables. Independent t-tests were used to assess differences between the variables. Kaplan–Meier curve and log-rank tests were used for survival data analysis, and univariate Cox regression analysis was used for survival factor analysis. Multivariate Cox regression analysis was used to identify independent prognostic factors. The accuracy of the prognostic model was assessed by time-dependent receiver operating characteristic (ROC) analysis. A P-value of <0.05 was considered statistically significant.

Results

Identification of Differentially Expressed Autophagy-Related Gene and Functional Enrichment Analysis

The mRNA levels of genes in gliomas (n = 698) and brain tissues (n = 5) in TCGA cohort were examined, and these values were compared using the Wilcoxon signed-rank test. This analysis revealed 15,001 DEGs in glioma tissues compared with normal brain tissues [FDR < 0.05, |log2 fold-change (FC)| > 1]. Then, we obtained 73 DEARGs in glioma tissues compared with normal brain tissues after taking the intersection of DEGs and ARGs (Figure 1A).

FIGURE 1

Figure 1. Identification of differentially expressed autophagy-related genes (DEARGs) in glioma and biological function analysis. (A) Volcano plot of DEARGs in tumor and normal samples of The Cancer Genome Atlas (TCGA) dataset; the green dots represent down-regulated genes, and red dots represent up-regulated genes. (B) The functional enrichment map of pathways. Each node indicates an enriched term and is colored by its cluster ID. (C) Enriched pathways of the genes positively correlated with DEARGs.

To explore the potential functional characteristics of the DEARGs, we annotated their functions using Metascape. The results showed that these genes were mainly involved in autophagy, animal autophagy, apoptosis, apoptotic signaling pathways, and the response to oxidative stress (Figures 1B,C).

Identification of Prognostic Differentially Expressed Autophagy-Related Genes

To identify the PDEARGs, we performed a univariate Cox regression analysis of each DEARG expression in TCGA cohort. A total of 62 DEARGs were significantly correlated with OS in glioma patients (P < 0.05).

Training Cohort to Identify Risk Genes for Inclusion in the Risk Model

Considering the effect of PDEARGs on patient prognosis, we further screened the PDEARGs by constructing a Cox regression hazard model, selecting the potential risk genes, and eliminating the overfitting genes in the model by LASSO regression; 23 genes (Figures 2A,B) were thus obtained and further analyzed by multivariate Cox proportional risk regression analysis. Finally, 15 risk genes (BRIC5, NFE2L2, GABARAP, IKBKE, BID, MAPK3, FKBP1B, MAPK8IP1, PRKCQ, CX3CL1, NPC1, HSP90AB1, DAPK2, SUPT20H, and PTEN) were identified, and the risk scores for each case were calculated using their expression levels and regression coefficients (Figure 2C). Among these genes, FKBP1B, BIRC5, NFE2L2, and IKBKE were considered high-risk genes (poor prognosis), while MAPK3, NPC1, DAPK2, MAPK8IP1, PRKCQ, CX3CL1, HSP90AB1, SUPT20H, GABARAP, PTEN, and BID were considered low-risk genes (protective factors) for patients’ OS. To explore the potential functional characteristics of the risk genes, we annotated their functions using Metascape. The results showed that these genes are mainly involved in animal autophagy, apoptotic signaling pathways, pathways in cancer, and the regulation of cellular response to stress (Figures 2D,E).

FIGURE 2

Figure 2. Construction of the prognostic risk model and biological function analysis. (A) Cross-validation for tuning parameter screening in the least absolute shrinkage and selection operator (LASSO) regression model. (B) LASSO coefficient profiles of the common genes. (C) Risk genes in the prognostic risk model. (D) The functional enrichment map of pathways of risk genes. Each node indicates an enriched term and is colored by its cluster ID. (E) Enriched pathways of the genes positively correlated with risk genes.

We then used clinical specimens from the Human Protein Atlas⁴ to analyze the expression of the proteins encoded by the risk genes (no data were found for DAPK2, PRKCQ, MAPK8IP1, and GABARAP). MAPK3 and CX3CL1 were strongly positive; NFE2L2, HSP90AB1, and SUPT20H were moderately positive; and FKBP1B, BIRC5, NPC1, IKBKE, BID, and PTEN were weakly positive in glioma tissue relative to their expression levels in normal tissue (Supplementary Figures S1A–K). Kaplan–Meier survival curves were also constructed to evaluate the relationship between the expression levels of the prognostic genes and OS, and the results showed that the low-expression groups of FKBP1B, BIRC5, NFE2L2, and IKBKE (P < 0.05) had better prognosis (Supplementary Figures S2B,C,E,G) and that the high-expression groups of MAPK3, NPC1, DAPK2, MAPK8IP1, PRKCQ, CX3CL1, HSP90AB1, SUPT20H, GABARAP, PTEN, and BID (P < 0.05) had better prognosis (Supplementary Figures S2A,D,F,H–O).

Construction and Verification of the Prognostic Risk Model

To assess the significance of the risk genes in the prognosis of patients with glioma in the risk models, we used estimated regression coefficients and the expression levels of the risk genes to calculate the risk scores for each patient. The calculation formula is as follows:

$$\begin{array}{c}\text{Risk} \text{score=}(\text{0}.\text{972364596}\times \text{MAPK3} \text{expression})\\ \text{+}(\text{0}.\text{299641268}\times \text{FKBP1B} \text{expression})\\ \text{+}(\text{-0}.\text{364718845}\times \text{MAPK8IP1} \text{expression})\\ \text{+}(\text{-0}.\text{271491791}\times \text{PRKCQ} \text{expression})\\ \text{+}(\text{0}.\text{314270977}\times \text{BIRC5} \text{expression})\\ \text{+}(\text{-0}.\text{172912776}\times \text{CX3CL1} \text{expression})\\ \text{+}(\text{0}.\text{664977132}\times \text{NPC1} \text{expression})\\ \text{+}(\text{-0}.\text{257268535}\times \text{HSP90AB1} \text{expression})\\ \text{+}(\text{0}.\text{81239816}\times \text{NFE2L2} \text{expression})\\ \text{+}(\text{0}.\text{271699944}\times \text{DAPK2} \text{expression})\\ \text{+}(\text{-0}.\text{289118873}\times \text{SUPT20H} \text{expression})\\ \text{+}(\text{-0}.\text{545908599}\times \text{GABARAP} \text{expression})\\ \text{+}(\text{-0}.\text{570898637}\times \text{PTEN} \text{expression})\\ \text{+}(\text{0}.\text{300429145}\times \text{IKBKE} \text{expression})\\ \text{+}(\text{-0}.\text{398268499}\times \text{BID} \text{expression}).\end{array}$$

In the CGGA cohort, the risk score was calculated using the same risk genes and regression coefficients for validation. The heatmaps show an overview of the correlation between ARG expression and clinical characteristics in TCGA and CGGA cohorts (Figure 3).

FIGURE 3

Figure 3. Heatmaps of autophagy-related genes based on the risk score value in The Cancer Genome Atlas (TCGA) cohort (A) and the Chinese Glioma Genome Atlas Chinese Glioma Genome Atlas (CGGA) cohort (B).

Patients in TCGA cohort were divided into a high-risk group (n = 332) and a low-risk group (n = 333), according to the median risk score. We used a Kaplan–Meier curve based on the log-rank test to identify the difference in prognosis between the high-risk group and the low-risk group. The prognosis of patients in the high-risk group was worse than that in the low-risk group (P < 0.05) (Figure 4A).

FIGURE 4

Figure 4. The prognostic value of the risk signature in The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA) cohorts. (A) Kaplan–Meier (K–M) survival analysis of the risk score in TCGA cohort. (B,C) Time-dependent receiver operating characteristic (ROC) curve analysis of the risk score and traditional factors in TCGA cohort. (D) Kaplan–Meier (K–M) survival analysis of the risk signature in the CGGA cohort. (E,F) Time-dependent ROC curve analysis of the risk score and traditional factors in the CGGA cohort.

We then used time-dependent ROC curves to confirm the accuracy of the model’s predictions for 1–10 years. The ROC (AUC) of the prediction model was 0.906 over 1 year, 0.921 over 2 years, 0.924 over 3 years, 0.89 over 4 years, 0.874 over 5 years, 0.882 over 6 years, 0.855 over 7 years, 0.809 over 8 years, 0.822 over 9 years, and 0.836 over 10 years in TCGA training cohort (Figure 4B), all of which were significantly higher than the other factors (Figure 4C), emphasizing the superior predictive value of the prognostic risk model. We then ranked the risk scores of the patients in TCGA cohort and analyzed their distribution (Figure 5A). The survival status of each patient in TCGA cohort is shown in Figure 5B.

FIGURE 5

Figure 5. Distribution of the risk score (A,D) survival status (B,E) and expression level of risk genes (C,F) in the Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA) cohort.

We also generated heatmap to describe the expression patterns of risk genes in both the high-risk and low-risk groups (Figure 5C). In patients with high-risk scores in TCGA cohort, the high-risk genes (BIRC5, FKBP1B, IKBKE, and NFE2L2) were up-regulated, while the protective genes (BID, CX3CL1, GABARAP, HSP90AB1, MAPK3, MAPK8IP1, NPC1, PRKCQ, PTEN, and SUPT20H) were down-regulated. In patients with low risk scores, these risk genes showed opposite patterns of expression.

To validate the accuracy of the prognostic risk model, the patients in the CGGA cohort were also divided into two groups based on the median risk score. Next, we used Kaplan–Meier survival analysis to identify the prognostic differences between the high-risk group and the low-risk group. The Kaplan–Meier survival curves differed significantly between the two risk groups in the CGGA cohort (P < 0.05) (Figure 4D), with the patients’ survival rates being higher in the low-risk group than in the high-risk group throughout the follow-up period. ROC analyses were also performed for the CGGA cohort over 1–10 years; the AUCs for each year were, in order, 0.698, 0.741, 0.717, 0.721, 0.718, 0.712, 0.677, 0.675, 0.659, and 0.626 (Figure 4E), which were higher than for the traditional factors (Figure 4F). The risk score distribution, survival status, and risk gene expression in the CGGA cohort are shown in Figures 5D,E. Similar to the results in TCGA cohort, the most protective gene levels were higher and the most risk gene levels lower in the low-risk group than in the high-risk group (Figure 5F). These results indicate that the prognostic risk model is capable of precisely predicting the prognosis of glioma patients.

In addition, when the patients were stratified according to different clinicopathologic parameters, the autophagy signature remained a significant prognostic factor in both TCGA cohort (Supplementary Figures S3A–J) and the CGGA cohort (Supplementary Figures S3K–T).

Independent Prognostic Value of the Prognostic Risk Model

Supplementary Table S1 shows the demographics and clinicopathologic characteristics of glioma patients in TCGA cohort and CGGA cohorts based on the autophagy signature. Univariate and multivariate Cox regression analyses were then performed to assess whether our model-generated risk scores were independent of other clinical parameters [WHO grade, age, gender, isocitrate dehydrogenase (IDH) mutation status, 1p/19q co-deletion status, and risk score] as prognostic factors for patients with glioma. Univariate analysis showed that WHO grade, age, IDH mutant status, 1p/19q co-deletion status, and risk score were associated with the prognosis of glioma patients in both TCGA and CGGA cohorts (P < 0.05) (Figures 6A,C); only gender was not. Multivariate analysis showed that age and risk score were independently associated with OS in TCGA cohort (P < 0.05), while WHO grade, 1p/19q co-deletion status, and risk score were independently associated with OS in the CGGA cohort (P < 0.05) (Figures 6B,D). These results suggest that the prognostic risk model can be independently used to predict the prognosis of patients with glioma.

FIGURE 6

Figure 6. Cox regression analyses in The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA) cohorts. (A,B) Univariate and multivariate Cox regression analyses of risk score and pathologic characteristics in TCGA cohort. (C,D) Univariate and multivariate Cox regression analyses of risk score and pathologic characteristics in the CGGA cohort.

Construction and Validation of the Nomogram

To better predict the prognosis of patients with glioma at 1–10 years, we constructed a new nomogram (Figure 7A). In Figures 7B,C, the blue lines represent the observed survival rate, the gray lines represent the ideal survival rate, and the black lines represent the optimized modified survival rate. The optimism-corrected line, also known as the bias-corrected or overfitting-corrected line, was produced using a bootstrap approach to estimate the observed and predicted values based on a non-parametric smoother applied to a sequence of predicted values (23). For the predicted 1- to 10-year survival plots in the training cohort (Figure 7B) and the validation cohort (Figure 7C), the observed and optimism-corrected lines are aligned, although the two are slightly different from the ideal 45° line, which means that the 1- to 10-year survival predicted by the nomogram reflects the observed actual survival in both the training and validation cohorts.

FIGURE 7

Figure 7. Nomogram for predicting overall survival of patients with glioma. (A) A nomogram integrating the signature risk score with the pathologic characteristics. (B,C) The calibration curve for the nomogram in The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) cohorts.

Clinical Utility of the Prognostic Risk Model

By analyzing somatic mutation data from TCGA dataset, we explored differences in genomic changes between the high-risk and low-risk groups, as some genomic changes have a negative impact on glioma survival (24). By comparing the frequency of mutation occurrences, we found that there was some significant mutation enrichment between the two groups – IDH1, TP53, ATRX, and CIC mutations were more abundant in the low-risk group (Figure 8A). By examining the association between risk scores and clinical features in both TCGA and CGGA cohorts, we found that the risk score was positively correlated with the grade of glioma (Figures 8C,F) and that it was significantly lower in the 1p/19q co-deletion and IDH mutation samples (Figures 8B,D,E,G).

FIGURE 8

Figure 8. Association between genomic alterations, pathologic characteristics, and the autophagy-related signature. (A) Differential somatic mutation analysis between high-risk group and low-risk group in The Cancer Genome Atlas (TCGA) cohort. (B–D) Distribution of the risk score in glioma patients stratified by 1p/19q codeletion status, WHO grade, and isocitrate dehydrogenase (IDH) mutation status in TCGA cohort. (E–G) Distribution of the risk score in glioma patients stratified by 1p/19q codeletion status, WHO grade, and IDH mutation status in the Chinese Glioma Genome Atlas (CGGA) cohort.

To examine the performance of our model in predicting the progression of glioma, we analyzed the relationships between the clinical variables in the two cohorts and the risk genes. We found that the expressions of BRIC5, FKBP1B, IKBKE, and NFE2L2 were significantly greater in the 1p/19q non-co-deletion samples, while the expressions of BID, CX3CL1, DAPK2, GABARAP, HSP90AB1, MAPK3, MAPK8IP1, and SUPT20H were significantly lower in those samples in both cohorts (Supplementary Figure S4); the expressions of BRIC5, FKBP1B, IKBKE, and NFE2L2 were positively correlated, and the expressions of BID, CX3CL1, DAPK2, MAPK3, MAPK8IP1, and PRKCQ negatively correlated with the grade of glioma in both cohorts (Supplementary Figure S5); and the expressions of BRIC5, FKBP1B, IKBKE, and NFE2L2 were significantly greater in IDH wild-type samples, while the expressions of BID, CX3CL1, GABARAP, HSP90AB1, MAPK3, MAPK8IP1, PRKCQ, and SUPT20H were significantly lower in those samples in both cohorts (Supplementary Figure S6) (P < 0.05). These results demonstrate that the dysregulation of autophagy-related risk gene expression is associated with the development of glioma.

To determine whether the model could reflect the status of the tumor immune microenvironment in glioma patients, we analyzed the correlation between the risk scores and immune cell infiltration by using Pearson correlation analysis in both TCGA and CGGA cohorts. We found that the risk scores were correlated positively with the levels of CD8+ T cells; gamma delta T cells; M0, M1, and M2 macrophages; and activated dendritic cells. Meanwhile, they were negatively correlated with the levels of naive B cells, naive CD4+ T cells, follicular helper T cells, activated NK cells, monocytes, and resting mast cells, suggesting that the immune-related gene signature could indicate the level of infiltrating immune cells in gliomas to a certain extent (Figures 9A,D). We also found that risk scores were correlated positively with the expressions of the immune checkpoint (B7-H3, CTLA4, LAG3, PD-1, PD-L1, PD-L2, and TIM-3) (Figures 9B,E) and inflammatory factors (HLA-A, HLA-B, and HLA-C) by using Pearson correlation analysis in TCGA and CGGA cohorts (Figures 9C,F), indicating that the risk scores are also positively associated with an elevated level of immune exhaustion and that the immune-related gene signature might promote activation of inflammatory responses in glioma patients. In conclusion, the immune-related gene signature could promote the malignant progression of glioma by shaping the tumor immune microenvironment.

FIGURE 9

Figure 9. Pearson correlation analysis between the risk score and immune cell infiltration, immune checkpoints, and human leukocyte antigen (HLA) complex expressions. (A–C) Pearson correlation analysis between the risk score and the level of immune cell infiltration, the expressions of immune checkpoints, and the expressions of HLA complex in The Cancer Genome Atlas (TCGA) cohort. (D–F) Pearson correlation analysis between the risk score and the level of immune cell infiltration, the expressions of immune checkpoints, and the expressions of HLA complex in the Chinese Glioma Genome Atlas (CGGA) cohort. “X” means P > 0.05.

Discussion

Glioma is a progressive disease that requires an effective prognostic indicator for diagnosis and treatment. There is increasing evidence that autophagy plays a key role in the occurrence, progression, and therapeutic resistance of glioma (4, 8, 25). And previous studies have shown that autophagy inhibits its development. In the KRAS-driven glioblastoma mouse model, autophagy blocked Atg7, Atg13, or Ulk1 by shRNA; inhibited the occurrence and growth of tumors; and extended the survival of mice (26). Chio et al. reported that treatment of temozolomide (TMZ)-resistant human glioma cells with honokiol can reduce cell activity and induce apoptosis, and Zhang et al. (27) showed that ARG MAPK8IP1 overexpression and SH3GLB1 knockdown inhibited glioma cell progression and improved TMZ sensitivity. There have also been many studies on targeted autophagy therapy for glioma in recent years; for example, 2-propylpentanoic acid (VPA) combined with TMZ has been found to have significant antitumor effects on TMZ-resistant malignant glioma cells by inducing autophagy (28). ARGs are therefore a promising therapeutic target and prognostic predictor in glioma.

In this study, we identified prognostic ARGs and used them to construct a reliable model for predicting the OS of patients with glioma. First, we analyzed the expression of 232 ARGs in glioma and found that the expression of 62 DEARGs was related to OS. These results suggest that ARGs are an important prognostic factor for glioma patients.

Next, we evaluated these PDEARGs in the stratification of patient outcomes. Through a combination of LASSO regression and Cox regression analysis, we identified 15 PDEARGs of interest and used them to construct a Cox regression risk model. We also analyzed and verified the reliability and stability of the model. The results showed that the model can accurately discriminate patients with different survival outcomes. Univariate and multivariate Cox regression analyses indicated that the prognostic risk model could independently predict the prognosis of patients with glioma. In addition, nomogram analysis indicated that the accuracy of the model in predicting glioma prognosis could be improved in combination with other clinical features.

Nomograms have been widely used in clinical work for visual presentation (29, 30). In this study, we created a nomogram containing WHO grade, age, gender, IDH mutation status, 1p/19q co-deletion status, and risk score. The calibration diagram based on the training and validation cohorts showed that the actual survival rates of patients with glioma were very close to the predicted survival rates, indicating that the prediction performance of the nomogram was excellent. The visual scoring system could help doctors and patients to make more personalized survival predictions, which could help them to choose better treatment options.

The relationship of risk genes with clinical variables (WHO grade, IDH mutation status, and 1p/19q codeletion status) was also analyzed. We found that risk genes were associated with the progression of gliomas. Our model therefore has high clinical application value in predicting the development of glioma.

A large number of studies have examined the molecular biomarkers of glioma using advances in large-scale genome sequencing technology (31–33). BIRC5 is a member of the family of genes that inhibit apoptosis, and, in glioma, BIRC5 expression is positively correlated with tumor grade (34, 35) – a recent study revealed that BIRC5 affects the tumorigenicity of glioma cells by regulating p53 protein (36). IKBKE is involved in the malignant transformation and development of tumors as an oncogene (37). Lu et al., have suggested that IKBKE plays a pivotal role in regulating cell proliferation and invasion and the epithelial–mesenchymal transition of malignant glioma cells in vitro and in vivo by impacting the Hippo pathway (38). The overexpression of IKBKE is closely related to the stage of the glioma, and the invasion and migration of IKBKE cells are reduced after IKBKE is silenced with synthetic siRNAs (39). The tumor-promoting activity of NFE2L2 has been attributed to its own function; its activity is generally enhanced in glioblastoma cell lines and tumors, and low NFE2L2 expression inhibits proliferation and self-renewal of glioma stem cells (40, 41). These results are consistent with our conclusions that the expressions of BIRC5, IKBKE, and NFE2L2 are positively correlated with the progression of glioma.

In addition to BRIC5, IKBKE, and NFE2L2, the 15 identified genes include many related to invasion and metastasis, such as PRKCQ, MAPK3, and PTEN. PRKCQ is a member of protein kinase C (PKC), and Couldwell et al., found a correlation between PKC activity and glioma cell proliferation, as glioma cells with a fast development speed have high PCK activity (42). MAPK3 is an important signal transduction molecule in the ERK/MAPK pathway, and previous studies have found that the enhanced functional activity of MAPK3 plays a critical role in the development and progression of gastric cancer (43, 44). Wang et al. have also investigated how MAPK3 regulates apoptosis and invasion in gliomas through targeting by mir-483-5p (45). PTEN dysregulation in mice leads to the development of multiple tumors (46, 47), and it has been reported that the up-regulation of PTEN expression is related to the growth inhibition of glioma cells in vitro (48).

Previous studies have demonstrated that immune infiltration plays a crucial role in determining therapeutic effects in glioma and the prognosis of patients with glioma (49, 50). The interaction between glioma cells and macrophages could promote the proliferation and invasion of tumor cells (51), and Weenink et al., found that the level of CD8+ T cells in glioma was positively correlated with poor prognosis (52). Therefore, we also explored the relationship between immune cell infiltration and risk scores in the prognostic risk model and found that risk scores were correlated, to some extent, with immune cell infiltration.

Previous studies have also demonstrated that the expressions of immune checkpoint and HLA complex are an important factor in determining the response to treatment and the prognosis of patients with glioma (53–55). Lemke et al. investigated the role of B7-H3 in glioblastoma and elucidated its mechanism in the most fatal type of glioblastoma (56), and Machulla et al. found that the expression of HLA was positively correlated with the occurrence of gliomas (57). We also found that the risk scores were positively correlated with the expressions of B7-H3, CTLA4, LAG3, PD-1, PD-L1, PD-L2, TIM-3, HLA-A, HLA-B, and HLA-C. These results confirm the reliability of the model in predicting the immune microenvironment of glioma.

There are two main limitations of the present study. Firstly, clinical information downloaded from TCGA and CGGA databases is limited and incomplete, and secondly, the identified mechanism of ARGs affecting glioma development requires further in vivo and in vitro experimental studies.

In summary, we used 15 ARGs to construct a prognostic risk model that accurately predicted the prognosis of glioma patients. The risk score generated by this model can be used as an independent prognostic indicator to discriminate patients with different survival outcomes. The model can also predict the immune microenvironment of glioma to a certain extent. However, further studies are needed using a prospective, large-scale, multicenter clinical cohort to validate the prognostic model.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/, http://www.cgga.org.cn/.

Author Contributions

CZ, KY, and RL conceived and designed the study, edited the manuscript, and acted as overall guarantors. YX and RL carried out the analysis procedure, contributed to the study design and interpretation of the results, and contributed to the drafting of the manuscript. XL, ND, DW, and LH contributed to the preparation of the dataset and drafts of the manuscript. All the authors reviewed the manuscript.

Funding

This work was supported by grants from the National Key R&D Program of China (2018YFA0109800).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.571189/full#supplementary-material

Footnotes

1. ^ https://tcga-data.nci.nih.gov/tcga/
2. ^ http://www.cgga.org.cn/
3. ^ http://www.autophagy.lu/
4. ^ http://www.proteinatlas.org

References

1. Ostrom QT, Cioffi G, Gittleman H, Patil N, Waite K, Kruchko C, et al. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2012-2016. Neurooncology. (2019) 21(Suppl. 5):v1–100. doi: 10.1093/neuonc/noz150

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Wang Y, Jiang T. Understanding high grade glioma: molecular mechanism, therapy and comprehensive management. Cancer Lett. (2013) 331:139–46. doi: 10.1016/j.canlet.2012.12.024

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Rock K, McArdle O, Forde P, Dunne M, Fitzpatrick D, O’Neill B, et al. A clinical review of treatment outcomes in glioblastoma multiforme–the validation in a non-trial population of the results of a randomised Phase III clinical trial: has a more radical approach improved survival? Br J Radiol. (2012) 85:e729–33. doi: 10.1259/bjr/83796755

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Kondo Y, Kanzawa T, Sawaya R, Kondo S. The role of autophagy in cancer development and response to therapy. Nat Rev Cancer. (2005) 5:726–34. doi: 10.1038/nrc1692

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Towers CG, Thorburn A. Therapeutic targeting of autophagy. EBioMedicine. (2016) 14:15–23. doi: 10.1016/j.ebiom.2016.10.034

PubMed Abstract | CrossRef Full Text | Google Scholar

6. White E. The role for autophagy in cancer. J Clin Invest. (2015) 125:42–6. doi: 10.1172/jci73941

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Cheng Y, Ren X, Hait WN, Yang JM. Therapeutic targeting of autophagy in disease: biology and pharmacology. Pharmacol Rev. (2013) 65:1162–97. doi: 10.1124/pr.112.007120

PubMed Abstract | CrossRef Full Text | Google Scholar

8. White E, Karp C, Strohecker AM, Guo Y, Mathew R. Role of autophagy in suppression of inflammation and cancer. Curr Opin Cell Biol. (2010) 22:212–7. doi: 10.1016/j.ceb.2009.12.008

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Guo JY, Chen HY, Mathew R, Fan J, Strohecker AM, Karsli-Uzunbas G, et al. Activated Ras requires autophagy to maintain oxidative metabolism and tumorigenesis. Genes Dev. (2011) 25:460–70. doi: 10.1101/gad.2016311

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Xie K, Liang C, Li Q, Yan C, Wang C, Gu Y, et al. Role of ATG10 expression quantitative trait loci in non-small cell lung cancer survival. Int J Cancer. (2016) 139:1564–73. doi: 10.1002/ijc.30205

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Guo JY, Karsli-Uzunbas G, Mathew R, Aisner SC, Kamphorst JJ, Strohecker AM, et al. Autophagy suppresses progression of K-ras-induced lung tumors to oncocytomas and maintains lipid homeostasis. Genes Dev. (2013) 27:1447–61. doi: 10.1101/gad.219642.113

CrossRef Full Text | Google Scholar

12. Cj P, Hv E, Vijayakurup V, Menon GR, Nair S, Gopala S. High LC3/Beclin expression correlates with poor survival in glioma: a definitive role for autophagy as evidenced by in vitro autophagic flux. Pathol Oncol Res POR. (2019) 25:137–48. doi: 10.1007/s12253-017-0310-7

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Giatromanolaki A, Sivridis E, Mitrakas A, Kalamida D, Zois CE, Haider S, et al. Autophagy and lysosomal related protein expression patterns in human glioblastoma. Cancer Biol Ther. (2014) 15:1468–78. doi: 10.4161/15384047.2014.955719

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zhang CB, Bao ZS, Wang HJ, Yan W, Liu YW, Li MY, et al. Correlation of IDH1/2 mutation with clinicopathologic factors and prognosis in anaplastic gliomas: a report of 203 patients from China. J Cancer Res Clin Oncol. (2014) 140:45–51. doi: 10.1007/s00432-013-1519-9

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Wang Z, Wang Z, Zhang C, Liu X, Li G, Liu S, et al. Genetic and clinical characterization of B7-H3 (CD276) expression and epigenetic regulation in diffuse brain glioma. Cancer Sci. (2018) 109:2697–705. doi: 10.1111/cas.13744

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Hu H, Mu Q, Bao Z, Chen Y, Liu Y, Chen J, et al. Mutational landscape of secondary glioblastoma guides MET-targeted trial in brain tumor. Cell. (2018) 175:1665–78.e18. doi: 10.1016/j.cell.2018.09.038

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Qian Z, Li Y, Fan X, Zhang C, Wang Y, Jiang T, et al. Prognostic value of a microRNA signature as a novel biomarker in patients with lower-grade gliomas. J Neuro Oncol. (2018) 137:127–37. doi: 10.1007/s11060-017-2704-5

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Tang G, Yin W. Development of an immune infiltration-related prognostic scoring system based on the genomic landscape analysis of glioblastoma multiforme. Front Oncol. (2020) 10:154. doi: 10.3389/fonc.2020.00154

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nature Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (Oxford, England). (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Linden A, Yarnold PR. Modeling time-to-event (survival) data using classification tree analysis. J Evaluation Clin Pract. (2017) 23:1299–308. doi: 10.1111/jep.12779

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. (2010) 33:1–22.

Google Scholar

23. Gittleman H, Sloan AE, Barnholtz-Sloan JS. An independently validated survival nomogram for lower grade glioma. Neuro Oncol. (2019) 22:665–74. doi: 10.1093/neuonc/noz191

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Aibaidula A, Chen AK, Shi Z, Li Y, Zhang R, Yang R, et al. Adult IDH wild-type lower-grade gliomas should be further stratified. Neuro Oncol. (2017) 19:1327–37. doi: 10.1093/neuonc/nox078

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Mathew R, Kongara S, Beaudoin B, Karp CM, Bray K, Degenhardt K, et al. Autophagy suppresses tumor progression by limiting chromosomal instability. Genes Dev. (2007) 21:1367–81. doi: 10.1101/gad.1545107

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Gammoh N, Fraser J, Puente C, Syred HM, Kang H, Ozawa T, et al. Suppression of autophagy impedes glioblastoma development and induces senescence. Autophagy. (2016) 12:1431–9. doi: 10.1080/15548627.2016.1190053

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Chio CC, Chen KY, Chang CK, Chuang JY, Liu CC, Liu SH, et al. Improved effects of honokiol on temozolomide-induced autophagy and apoptosis of drug-sensitive and -tolerant glioma cells. BMC Cancer. (2018) 18:379. doi: 10.1186/s12885-018-4267-z

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Zhang H, Lu X, Wang N, Wang J, Cao Y, Wang T, et al. Autophagy-related gene expression is an independent prognostic indicator of glioma. Oncotarget. (2017) 8:60987–1000. doi: 10.18632/oncotarget.17719

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Harrell FE Jr., Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Stat Med. (1996) 15:361–87. doi: 10.1002/(sici)1097-0258(19960229)15:43.0.Co;2-4

CrossRef Full Text | Google Scholar

30. Liu Y, Wu L, Ao H, Zhao M, Leng X, Liu M, et al. Prognostic implications of autophagy-associated gene signatures in non-small cell lung cancer. Aging. (2019) 11:11440–62. doi: 10.18632/aging.102544

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Killela PJ, Pirozzi CJ, Healy P, Reitman ZJ, Lipp E, Rasheed BA, et al. Mutations in IDH1, IDH2, and in the TERT promoter define clinically distinct subgroups of adult malignant gliomas. Oncotarget. (2014) 5:1515–25. doi: 10.18632/oncotarget.1765

CrossRef Full Text | Google Scholar

32. Sturm D, Witt H, Hovestadt V, Khuong-Quang DA, Jones DT, Konermann C, et al. Hotspot mutations in H3F3A and IDH1 define distinct epigenetic and biological subgroups of glioblastoma. Cancer Cell. (2012) 22:425–37. doi: 10.1016/j.ccr.2012.08.024

CrossRef Full Text | Google Scholar

33. Chai RC, Zhang KN, Chang YZ, Wu F, Liu YQ, Zhao Z, et al. Systematically characterize the clinical and biological significances of 1p19q genes in 1p/19q non-codeletion glioma. Carcinogenesis. (2019) 40:1229–39. doi: 10.1093/carcin/bgz102

CrossRef Full Text | Google Scholar

34. Uematsu M, Ohsawa I, Aokage T, Nishimaki K, Matsumoto K, Takahashi H, et al. Prognostic significance of the immunohistochemical index of survivin in glioma: a comparative study with the MIB-1 index. J Neuro Oncol. (2005) 72:231–8. doi: 10.1007/s11060-004-2353-3

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Liu X, Chen N, Wang X, He Y, Chen X, Huang Y, et al. Apoptosis and proliferation markers in diffusely infiltrating astrocytomas: profiling of 17 molecules. J Neuropathol Exp Neurol. (2006) 65:905–13. doi: 10.1097/01.jnen.0000235857.79502.c3

CrossRef Full Text | Google Scholar

36. Conde M, Michen S, Wiedemuth R, Klink B, Schrock E, Schackert G, et al. Chromosomal instability induced by increased BIRC5/Survivin levels affects tumorigenicity of glioma cells. BMC Cancer. (2017) 17:889. doi: 10.1186/s12885-017-3932-y

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Boehm JS, Zhao JJ, Yao J, Kim SY, Firestein R, Dunn IF, et al. Integrative genomic approaches identify IKBKE as a breast cancer oncogene. Cell. (2007) 129:1065–79. doi: 10.1016/j.cell.2007.03.052

CrossRef Full Text | Google Scholar

38. Lu J, Yang Y, Guo G, Liu Y, Zhang Z, Dong S, et al. IKBKE regulates cell proliferation and epithelial-mesenchymal transition of human malignant glioma via the Hippo pathway. Oncotarget. (2017) 8:49502–14. doi: 10.18632/oncotarget.17738

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Li H, Chen L, Zhang A, Wang G, Han L, Yu K, et al. Silencing of IKKepsilon using siRNA inhibits proliferation and invasion of glioma cells in vitro and in vivo. Int J Oncol. (2012) 41:169–78. doi: 10.3892/ijo.2012.1452

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Cong ZX, Wang HD, Wang JW, Zhou Y, Pan H, Zhang DD, et al. ERK and PI3K signaling cascades induce Nrf2 activation and regulate cell viability partly through Nrf2 in human glioblastoma cells. Oncol Rep. (2013) 30:715–22. doi: 10.3892/or.2013.2485

CrossRef Full Text | Google Scholar

41. Kanamori M, Higa T, Sonoda Y, Murakami S, Dodo M, Kitamura H, et al. Activation of the NRF2 pathway and its impact on the prognosis of anaplastic glioma patients. Neuro Oncol. (2015) 17:555–65. doi: 10.1093/neuonc/nou282

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Couldwell WT, Antel JP, Yong VW. Protein kinase C activity correlates with the growth rate of malignant gliomas: Part II. Effects of glioma mitogens and modulators of protein kinase C. Neurosurgery. (1992) 31:717–24. doi: 10.1227/00006123-199210000-00015

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Wang X, Yu Z, Zhou Q, Wu X, Chen X, Li J, et al. Tissue transglutaminase-2 promotes gastric cancer progression via the ERK1/2 pathway. Oncotarget. (2016) 7:7066–79. doi: 10.18632/oncotarget.6883

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Li PY, Lv J, Qi WW, Zhao SF, Sun LB, Liu N, et al. Tspan9 inhibits the proliferation, migration and invasion of human gastric cancer SGC7901 cells via the ERK1/2 pathway. Oncol Rep. (2016) 36:448–54. doi: 10.3892/or.2016.4805

CrossRef Full Text | Google Scholar

45. Wang L, Shi M, Hou S, Ding B, Liu L, Ji X, et al. MiR-483-5p suppresses the proliferation of glioma cells via directly targeting ERK1. FEBS Lett. (2012) 586:1312–7. doi: 10.1016/j.febslet.2012.03.035

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Shang Y, Guo XX, Li WW, Rao W, Chen ML, Mu LN, et al. Cucurbitacin-B inhibits neuroblastoma cell proliferation through up-regulation of PTEN. Eur Rev Med Pharmacol Sci. (2014) 18:3297–303.

Google Scholar

47. Li XT, Wang HZ, Wu ZW, Yang TQ, Zhao ZH, Chen GL, et al. miR-494-3p regulates cellular proliferation, invasion, migration, and apoptosis by PTEN/AKT signaling in human glioblastoma cells. Cell Mol Neurobiol. (2015) 35:679–87. doi: 10.1007/s10571-015-0163-0

CrossRef Full Text | Google Scholar

48. Tu M, Cai L, Zheng W, Su Z, Chen Y, Qi S. CD164 regulates proliferation and apoptosis by targeting PTEN in human glioma. Mol Med Rep. (2017) 15:1713–21. doi: 10.3892/mmr.2017.6204

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Pereira MB, Barros LRC, Bracco PA, Vigo A, Boroni M, Bonamino MH, et al. Transcriptional characterization of immunological infiltrates and their relation with glioblastoma patients overall survival. Oncoimmunology. (2018) 7:e1431083. doi: 10.1080/2162402x.2018.1431083

CrossRef Full Text | Google Scholar

50. Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med. (2015) 21:938–45. doi: 10.1038/nm.3909

CrossRef Full Text | Google Scholar

51. Morisse MC, Jouannet S, Dominguez-Villar M, Sanson M, Idbaih A. Interactions between tumor-associated macrophages and tumor cells in glioblastoma: unraveling promising targeted therapies. Expert Rev Neurother. (2018) 18:729–37. doi: 10.1080/14737175.2018.1510321

CrossRef Full Text | Google Scholar

52. Weenink B, Draaisma K, Ooi HZ, Kros JM, Smit PAES, Debets R, et al. Low-grade glioma harbors few CD8 T cells, which is accompanied by decreased expression of chemo-attractants, not immunogenic antigens. Sci Rep. (2019) 9:14643. doi: 10.1038/s41598-019-51063-6

CrossRef Full Text | Google Scholar

53. Belcaid Z, Phallen JA, Zeng J, See AP, Mathios D, Gottschalk C, et al. Focal radiation therapy combined with 4-1BB activation and CTLA-4 blockade yields long-term survival and a protective antigen-specific memory response in a murine glioma model. PLoS One. (2014) 9:e101764. doi: 10.1371/journal.pone.0101764

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Cloughesy TF, Mochizuki AY, Orpilla JR, Hugo W, Lee AH, Davidson TB, et al. Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nat Med. (2019) 25:477–86. doi: 10.1038/s41591-018-0337-7

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Feng E, Liang T, Wang X, Du J, Tang K, Wang X, et al. Correlation of alteration of HLA-F expression and clinical characterization in 593 brain glioma samples. J Neuroinflamm. (2019) 16:33. doi: 10.1186/s12974-019-1418-3

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Lemke D, Pfenning PN, Sahm F, Klein AC, Kempf T, Warnken U, et al. Costimulatory protein 4IgB7H3 drives the malignant phenotype of glioblastoma by mediating immune escape and invasiveness. Clin Cancer Res. (2012) 18:105–17. doi: 10.1158/1078-0432.Ccr-11-0880

CrossRef Full Text | Google Scholar

57. Machulla HK, Steinborn F, Schaaf A, Heidecke V, Rainov NG. Brain glioma and human leukocyte antigens (HLA)–is there an association. J Neuro Oncol. (2001) 52:253–61. doi: 10.1023/a:1010612327647

CrossRef Full Text | Google Scholar