An immune signature to predict the prognosis of ATRX-wildtype glioma patients and guide immune checkpoint blockade therapy

High- and low-immune/stromal grouping of ATRX-mt and ATRX-wt glioma patients in TCGA

Our workflow to identify, test and validate prognostic and predictive biomarkers of glioma is shown in Figure 1A. The ESTIMATE algorithm was used to calculate immune and stromal scores for ATRX-mutation-type (ATRX-mt) and ATRX-wt glioma patients from TCGA, and the median scores were used to divide patients into high- and low-scoring groups. According to a Kaplan-Meier survival analysis, high vs. low immune/stromal scores did not significantly impact ATRX-mt glioma patient survival (Supplementary Table 1; Figure 1B, 1C). However, high immune and stromal scores were associated with poorer overall survival than low immune/stromal scores in ATRX-wt glioma patients (Figure 1D, 1E).

Figure 1. Stromal and immune scores were associated with survival in ATRX-wt glioma patients. (A) Flow chart of the analytical process in this study. (B) Kaplan-Meier survival analysis for ATRX-mt glioma patients in the high- and low-immune-score groups. (C) Kaplan-Meier survival analysis for ATRX-mt glioma patients in the high- and low-stromal-score groups. (D) Kaplan-Meier survival analysis for ATRX-wt glioma patients in the high- and low-immune-score groups. (E) Kaplan-Meier survival analysis for ATRX-wt glioma patients in the high- and low-stromal-score groups.

Next, we compared the gene expression profiles of ATRX-wt glioma patients with high and low immune/stromal scores. We found 166 upregulated and 28 downregulated genes in the high-stromal-score group compared with the low-stromal-score group (Figure 2A, 2B). In addition, we identified 158 upregulated and 44 downregulated genes in high-immune-score patients compared with low-immune-score patients (Figure 2C, 2D). When we analyzed the differentially expressed genes (DEGs) that overlapped between the high-stromal-score and high-immune-score groups of ATRX-wt glioma tissues, we found 136 upregulated and 26 downregulated genes (Supplementary Table 2; Figure 2E, 2F).

Figure 2. DEGs in ATRX-wt glioma tissues with high vs. low stromal/immune scores. (A) Volcano plot showing the DEGs between the high- and low-stromal-score groups of ATRX-wt glioma tissues. (B) Heat map showing the DEGs between the high- and low-stromal-score groups of ATRX-wt glioma tissues. (C) Volcano plot showing the DEGs between the high- and low-immune-score groups of ATRX-wt glioma tissues. (D) Heat map showing the DEGs between the high- and low-immune-score groups of ATRX-wt glioma tissues. (E) Overlapping downregulated genes between the high-stromal-score and high-immune-score groups of ATRX-wt glioma tissues. (F) Overlapping upregulated genes between the high-stromal-score and high-immune-score groups of ATRX-wt glioma tissues.

We then generated a protein-protein interaction (PPI) network based on the overlapping DEGs in ATRX-wt glioma patients with high stromal and immune scores. The PPI network contained 80 interacting genes (Figure 3A). Candidate hub genes from this network were subjected to Gene Ontology analyses for Molecular Function, Biological Process and Cellular Component. The Molecular Function terms associated with the candidate hub genes were enriched in “RNA polymerase II transcription factor activity”, “RNA polymerase II core promoter proximal region sequence specific”, “sequence specific double stranded”, “DNA binding” and “RNA polymerase II transcription regulatory region sequence specific binding” (Figure 3B). In the Biological Process analysis, the candidate hub genes were found in “regulation of transcription from RNA polymerase II promoter”, “positive regulation of transcription from RNA polymerase II promoter”, “anterior/posterior pattern specification”, “immune response” and “embryonic skeletal system morphogenesis” (Figure 3C). The Cellular Component analysis was enriched in “cell surface”, “external side of plasma membrane”, “nucleoplasm”, “chromatin” and “nucleus” (Figure 3D). We also performed a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, which revealed that the candidate hub genes were involved in the “T cell receptor signaling pathway”, “chemokine signaling pathway”, “intestinal immune network for IgA production”, “cell adhesion molecules” and “IL-17 signaling pathway” (Figure 3E).

Figure 3. Landscape of the 162 overlapping DEGs. (A) A PPI network was constructed using the 162 overlapping DEGs, and isolated genes were removed. Genes in the PPI network were set as candidate hub genes. (B) Molecular Function analysis of candidate hub genes. (C) Biological Process analysis of candidate hub genes. (D) Cellular Component analysis of candidate hub genes. (E) KEGG analysis of candidate central genes.

Constructing the immune profiles of ATRX-wt glioma patients

Next, we performed a univariate Cox regression analysis, which revealed that seven of the 80 aforementioned candidate genes were associated with survival in ATRX-wt glioma patients (Table 1). We subsequently conducted a LASSO analysis (Figure 4A, 4B) and a multivariate Cox analysis (Figure 4C), which identified homeobox A5 (HOXA5), protein tyrosine phosphatase non-receptor type 2 (PTPN2I, WT1 transcription factor (WT1), homeobox D10 (HOXD10), periostin (POSTN), ADAM like decysin 1 (ADAMDEC1) and myosin binding protein H (MYBPH) as independent predictors of survival in ATRX-wt glioma patients. These results were used to develop a risk model capable of predicting survival in ATRX-wt glioma patients. After computer optimization, and the following risk score expression was obtained: 0.11 × HOXA5 expression + 0.41 × PTPN2 expression + 0.14 × WT1 expression + 0.08 HOXD10 expression + 0.14 × POSTN expression + 0.11 × ADAMDEC1 expression + 0.10 × MYBPH expression.

Figure 4. Key genes selected for risk model construction. (A, B) LASSO analysis of key genes associated with survival in ATRX-wt glioma patients. (C) Multivariate Cox regression analysis of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH. These seven genes were used for risk model construction.

Risk model validation for ATRX-wt glioma patients from TCGA

Next, we randomized ATRX-wt glioma patients from TCGA into training and testing groups to evaluate the applicability of the risk model. The median risk score of ATRX-wt glioma patients in the training cohort was used to divide patients into high-risk and low-risk groups (Figure 5A). The overall survival rate was lower in the high-risk group than in the low-risk group (Figure 5B). The area under the curve (AUC) values for predicting survival after one, three and five years in ATRX-wt glioma patients from TCGA were 0.905, 0.917 and 0.883, respectively (Figure 5C–5E). The proportion of deaths among ATRX-wt glioma patients in the training cohort was higher in the high-risk-score group than in the low-risk-score group (Figure 5F). In glioma tissues from the high-risk group of the training cohort, HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH were highly expressed (Figure 5G).

Figure 5. Validation of the applicability of the risk model in the training cohort of ATRX-wt glioma patients. (A) ATRX-wt glioma patients in the training cohort were divided into high- and low-risk-score groups based on the median risk score. (B) Survival differences between ATRX-wt glioma patients in the high- and low-risk-score groups in the training cohort. (C–E) Prognostic value of the risk model for the one-, three- and five-year survival of ATRX-wt glioma patients in the training cohort. (F) Deaths of ATRX-wt glioma patients in the high- and low-risk-score groups in the training cohort (green dots represent living cases; red dots represent dead cases). (G) Expression of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH in glioma patients in the high- and low-risk-score groups in the training cohort.

We then validated the model using the test cohort. Again, ATRX-wt glioma patients were divided into high- and low-risk groups based on the median risk score (Figure 6A). Glioma patients with high-risk scores were more likely to die than those with low-risk scores (Figure 6B). The AUCs for predicting survival after one, three and five years in ATRX-wt glioma patients from the test cohort were 0.882, 0.885 and 0.825, respectively (Figure 6C–6E). The death rate of the high-risk group of glioma patients was higher in the test cohort than in the high-risk-score group than in the low-risk-score group (Figure 6F). In the test cohort, HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH levels were elevated in glioma tissues with high-risk scores (Figure 6G). Thus, our risk model exhibited significant prognostic value in ATRX-wt glioma patients.

Figure 6. Validation of the applicability of the risk model in the test cohort of ATRX-wt glioma patients. (A) ATRX-wt glioma patients in the test cohort were divided into high- and low-risk-score groups based on the median risk score. (B) Survival differences between ATRX-wt glioma patients in the high- and low-risk-score groups in the test cohort. (C–E) Prognostic value of the risk model for one-, three- and five-year survival in ATRX-wt glioma patients in the test cohort. (F) Deaths of ATRX-wt glioma patients in the high- and low-score groups in the test cohort (green dots represent living cases; red dots represent dead cases). (G) Expression of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH in glioma patients in the high- and low-risk-score groups of the test cohort.

Exploring the applicability of the risk model for ATRX-mt glioma patients from TCGA

Next, we classified ATRX-mt patients from TCGA into high- and low-risk groups based on the median risk score (Figure 7A). High risk scores were associated with a lower survival rate than low risk scores (Figure 7B). However, the AUCs for predicting one-, three- and five-year survival among ATRX-mt glioma patients were 0.53, 0.543 and 0.524, respectively (Figure 7C–7E). Moreover, the death rates for ATRX-mt glioma patients did not differ significantly between the high- and low-risk groups (Figure 7F). Nevertheless, HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH were overexpressed in ATRX-mt glioma tissues with high-risk scores compared with low-risk scores (Figure 7G). These results demonstrated that the risk model constructed with HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH could not accurately predict survival in ATRX-mt glioma patients, although it could in ATRX-wt glioma patients.

Figure 7. Validation of the applicability of the risk model in ATRX-mt glioma patients from TCGA. (A) ATRX-mt glioma patients from TCGA were divided into high- and low-risk-score groups based on the median risk score. (B) Survival differences between the high- and low-risk-score groups of ATRX-mt glioma patients from TCGA. (C–E) Prognostic value of the risk model for the one-, three- and five-year survival of ATRX-mt glioma patients from TCGA. (F) Deaths of ATRX-mt glioma patients from TCGA in the high- and low-risk-score groups (green dots represent living cases; red dots represent dead cases). (G) Expression of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH in ATRX-mt glioma tissues from the high- and low-risk-score groups.

Immune characteristics as an independent prognostic factor for ATRX-wt glioma patients

Our multivariate Cox regression analysis indicated that the HOXA5-derived immune signature was independently associated with the outcomes of ATRX-wt glioma patients. We used this signature risk score and patients’ clinical characteristics to create a nomogram (Figure 8A), which had high prognostic value for survival at one, three and five years (Figure 8B).

Figure 8. Construction of the nomogram. (A) Construction of a nomogram using age, sex, grade and risk score. (B) Nomogram in glioma patients with one-, three- and five-year survival.

Expression of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH in ATRX-wt glioma tissues

Next, we obtained 54 glioma tissues from patients at Guizhou Medical University Affiliated Hospital, and divided them according to whether the patients survived long-term (≥15 months) or short-term (<15 months). Then, we performed quantitative real-time PCR and immunohistochemical analyses to detect the mRNA and protein levels of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH in these tissues. These genes were expressed at significantly higher levels in glioma tissues from short-term survivors than in those from long-term survivors (Supplementary Figure 1; Figure 9A, 9B). A receiver operating characteristic curve analysis indicated that HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH had significant prognostic value for the survival of ATRX-wt glioma patients (AUC = 0.84, 0.81, 0.79, 0.91, 0.80, 0.79 and 0.85, respectively) (Figure 9C).

Figure 9. Expression of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH in ATRX-wt glioma tissues. ATRX-wt patients were divided into long-term and short-term survival groups based on a cut-off of 15 months. (A) Immunohistochemistry scores for HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH in ATRX-wt glioma tissues. (B) Representative plots of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH expression in ATRX-wt glioma tissues from the long- and short-term survival groups. (C) Prognostic value of HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH1 in distinguishing ATRX-wt glioma patients with long- and short-term survival.

Immunological properties of immune markers

We then used the “CIBERSORT” R package to evaluate immune cell infiltration in ATRX-wt glioma tissues from TCGA, and found that 22 types of immune cells infiltrated differently in the high- vs. low-risk groups. CD8+ T cells, naive CD4+ T cells, monocytes, M1 macrophages and M0 macrophages were in higher proportion in glioma tissues with low-risk scores (Figure 10A–10C). We also evaluated the response to ICB treatment in the high- and low-risk ATRX-wt groups, and found that patients in the low-risk group had a higher response rate to ICB treatment (Figure 10D). Thus, our risk model can be used to predict the prognosis and treatment responsiveness of ATRX-wt patients.

Figure 10. Immunological characteristics of the three immune features. (A, B) Gene expression profiles of the high- and low-risk-score groups of ATRX-wt glioma tissues from TCGA were transformed into 22 immune cell expression matrices. (C) Immune cell differences between the high- and low-risk groups of ATRX-wt glioma tissues from TCGA. (D) Responders and non-responders to ICB treatment among ATRX-wt glioma patients in the high- and low-risk groups. ^*P < 0.05; ^**P < 0.01; ^***P < 0.001.

Patients with high-risk glioma may benefit from rapamycin, dasatinib, 5-fluorouracil and gemcitabine

Lastly, we used the OncoPredict algorithm to derive a drug sensitivity model from the gene expression data in the ATRX-wt atlas. We evaluated a total of 198 inhibitors, and noted that high-risk glioma tissues were predicted to be sensitive to rapamycin, dasatinib, 5-fluorouracil and gemcitabine (Supplementary Table 3; Figure 11A–11E). This evidence suggests that rapamycin, dasatinib, 5-fluorouracil and gemcitabine may be useful drugs for patients with high-risk gliomas.

Figure 11. Selection of appropriate drugs for glioma patients in the high-risk group. (A–E) OncoPredict showed that the drug scores for rapamycin, dasatinib, 5-fluorouracil and gemcitabine differed between glioma patients in the high- and low-risk groups. ^***P < 0.001.

Downloading and preprocessing gene expression profiles

The gene expression profiles and clinical characteristics of 452 ATRX-wt and 188 ATRX-mt glioma patients were obtained from TCGA. The gene expression profiles were normalized and centralized, and gene names were assigned to the probes. The immune and stromal scores of the ATRX-wt glioma tissues were calculated using the ESTIMATE algorithm.

DEG analysis

The median immune and stromal scores from the dataset in TCGA were used to divide patients into high- or low-scoring groups. DEG analysis was performed using EdgeR, with an adjusted P-value of 0.05 and a |logFold-Change| <1 set as the threshold for significance. Volcano plots and heat maps were used to visualize gene expression changes in the groups with high immune/stromal scores [34].

PPI network

A PPI network was constructed by mapping DEG information to the Search Tool for the Retrieval of Interacting Genes (STRING) database. Isolated genes were removed from the original PPI network using Cytoscape software. A visual analysis was performed, and reciprocally related genes were designated as hub genes and included in the next step.

Functional and pathway enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery was used to analyze the enriched KEGG and Gene Ontology terms of genes. Three categories were used for the Gene Ontology analysis: Biological Processes, Cellular Components and Molecular Functions. A bubble diagram was used to present the terms to determine significance.

Immune signature construction and verification

In ATRX-wt patients, a univariate Cox regression analysis with a significance threshold was used to construct an immune signature. A LASSO operator with an appropriate penalty was used to eliminate genes with the same genetic information. A prognostic risk model was developed using multivariate Cox regression analysis with the Akaike information criterion. Kaplan-Meier survival analyses and receiver operating characteristic curve analyses were used for ATRX-wt and ATRX-mt glioma patients [35].

Construction of column line diagrams

The “rms” package in R was used to create column line plots and conduct a multi-factor regression analysis. Asymptotic lines were then used to plot on the same plane at a certain scale. The accuracy of the line plot was predicted, and the prognostic value of the line plot was determined.

Immune cell analysis

Using the “CIBERSORT” R package, we examined 22 immune cells infiltrating ATRX-wt glioma tissues in TCGA. Differentially infiltrated cells in the high- and low-risk-score groups were analyzed using the unpaired t-test, with significance set as P < 0.05.

Immunohistochemical analysis

Throughout Guizhou Medical University Affiliated Second People Hospital of Guiyang, 54 ATRX-wt glioma tissues were collected with approval from the Guizhou Medical University Human Ethics Committee. All participants provided informed consent before they were given radiotherapy or chemotherapy. Immunohistochemical staining was performed as described in a previous study [36]. The sections were probed with the following antibodies: HOXA5 (1:200; ab140636, Abcam, Cambridge, UK), PTPN2 (1:100; 11214-1-AP, Proteintech, Wuhan, China), WT1 (1:100; 12609-1-AP, Proteintech), POSTN (1:100; 66491-1-Ig, Proteintech), ADAMDEC1 (1:100; 17899-1-AP, Proteintech), MYBPH (1:100; ab197216, Abcam) and HOXD10 (1:100; ab138508, Abcam).

Drug score analysis

In vivo drug responses can be predicted using OncoPredict, an algorithm developed by Maeser et al. [37]. In order to calculate the drug sensitivity of gliomas, OncoPredict scripts were used to match the gene expression matrix of each glioma sample to the chemotherapeutic effects of drugs recorded in Cancer and the gene expression information for cancer lines in the Broad Institute Cancer Cell Line Encyclopedia. Glioma patients with higher drug scores are less sensitive to drugs. The limma package was used to analyze the differences in drug scores between those at high and low risk. A |logFold-Change| ≥1 and adjusted P < 0.05 were used as cut-offs for determining significance.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.