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

Immune and stromal cells contribute to glioma progression by infiltrating the tumor microenvironment. We used clinical characteristics, RNA sequencing data and the ESTIMATE algorithm to obtain stromal and immune scores for alpha thalassemia retardation syndrome X-linked (ATRX)-mutation-type (ATRX-mt) and ATRX-wildtype (ATRX-wt) glioma tissues from The Cancer Genome Atlas. To identify specific immune biomarkers of glioma, we compared the gene expression profiles of ATRX-wt glioma tissues with high vs. low immune/stromal scores, and discovered 162 differentially expressed genes. The protein-protein interaction network based on these results contained 80 interacting genes, of which seven (HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH) were identified as key prognostic genes via LASSO and Cox regression analyses. A risk model constructed using the expression of these seven genes could predict survival for ATRX-wt glioma patients, but was ineffective for ATRX-mt patients. T cells and macrophages were more prevalent in low-risk than in high-risk glioma tissues. Immune checkpoint blockade treatment was highly beneficial for patients with low risk scores. High-risk gliomas were predicted to be more sensitive to rapamycin, dasatinib, 5-fluorouracil and gemcitabine. Thus, our model can be used for the diagnosis, prognostic prediction and treatment planning of ATRX-wt glioma patients.


INTRODUCTION
Among primary malignant tumors of the nervous system, gliomas are the most common.In recent years, molecular technology has revealed some of the genetic and chromosomal changes associated with the occurrence, development and prognosis of gliomas [1].Gliomas can be characterized by several pathological molecular changes, including telomerase reverse transcriptase promoter mutations, isocitrate dehydrogenase mutations, methylguanine methyltransferase promoter methylation, 1p/19q codeletion and alpha thalassemia retardation syndrome X-linked (ATRX) mutations [2].AGING At tandem repeat sequences in the genome, ATRX helps to prevent replication fork arrest, facilitate histone variant formation and prevent homologous recombination at telomeres [3].Loss of ATRX causes epigenetic changes, including hypomethylation of repetitive elements such as telomeres [4].ATRX mutation or loss is common in a variety of tumor types, including low-grade astrocytomas and secondary glioblastomas [5].A bioinformatic analysis revealed that glioblastoma patients with ATRX loss had longer overall survival and greater benefit from temozolomide treatment than patients without this change [6].However, it is unclear how the dysregulation of stromal and immune cell infiltration in ATRX-wildtype (ATRX-wt) glioma influences its development and progression.Therefore, identifying specific biomarkers of ATRX-wt glioma may facilitate the treatment of this disease.
The Cancer Gene Atlas (TCGA) is commonly used to identify tumor biomarkers in bioinformatic analyses [7].Zhu et al. found that a nuclear translocation-associated gene signature combined with isocitrate dehydrogenase mutation status and 1p/19q codeletion status could improve prognostic prediction in glioma patients [8].Feng et al. identified and validated an autophagy-associated signature to predict survival in low-grade gliomas [9].In this study, we used TCGA data to construct a specific risk model that predicts the prognosis of ATRX-wt glioma patients and informs immune checkpoint blockade (ICB) therapy.Our results may provide new insights into the diagnosis and treatment of ATRX-wt gliomas.

RESULTS
High-and low-immune/stromal grouping of ATRXmt 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-mutationtype (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 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-stromalscore and high-immune-score groups of ATRX-wt glioma tissues, we found 136 upregulated and 26 downregulated genes (Supplementary Table 2; Figure 2E,  2F).
We then generated a protein-protein interaction (PPI) network based on the overlapping DEGs in ATRXwt 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).

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

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

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

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

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

DISCUSSION
Abnormal ATRX expression has been detected in a variety of malignant tumors [10,11].ATRX mutations have been observed in neuroblastomas, pancreatic neuroendocrine tumors and pediatric osteosarcomas [12,13].ATRX protein contains a C-terminal helicase/ATPase domain, and belongs to the SWI/SNF2 chromatin remodeling protein family [14].The N-terminal ATRX-DNA methyltransferase 3-DNA methyltransferase 3-like ('ADD') domain contains a plant homeodomain and a GATA zinc finger structure [15].The GATA zinc finger can bind to DNA or chromatin, while the plant homeodomain participates in chromatin regulation and transcription [16,17].
Previous studies have indicated that the immune microenvironment influences the progression of ATRXwt gliomas [12].We compared the survival of ATRXwt glioma patients with high and low stromal/immune scores, and found that patients with low scores had a higher survival rate.We also assessed the DEGs between ATRX-wt glioma patients with high and low stromal/ immune scores, and used them to generate a PPI network.We found that 80 of the 162 DEGs interacted with other genes and were significantly associated with patients' prognoses.Using LASSO and Cox regression analyses, we identified HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH as independent predictors of survival.Based on the expression of these genes, we then generated immune profiles to classify ATRX-wt glioma patients as high-or low-risk.
The seven genes identified in this study have been associated with cancer in previous studies.For instance, alterations in the HOX family members HOXA5 and HOXD10 have been implicated in the development and progression of cancer [18,19].In non-small cell lung cancer, HOXA5 was found to promote apoptosis and inhibit proliferation by upregulating linc00312 expression [20].HOXD10 was identified as a biological correlate of tumor suppressor DNA and an inducer of miRNA-7 and insulin-like growth factor binding protein 3 in colorectal cancer [21].In tumor cells, PTPN2 was shown to enhance antigen presentation and growth arrest [22].WT1 has been detected in hematologic malignancies and solid tumors (breast, lung, pancreatic and prostate cancers).Furthermore, WT1 protein has high immunogenicity, suggesting that it may be a useful therapeutic agent in patients with WT1 gene amplification [23].Ovarian cancer cells incorporating POSTN from cancer-associated fibroblasts were   reported to migrate and invade more effectively due to phosphoinositide 3-kinase/Akt pathway activation [24].In the same study, the pro-metastatic and fibroblastactivating properties of transforming growth factor β1 were shown to depend partially on POSTN [24].MYBPH was found to suppress Rho-associated coiledcoil containing protein kinase 1 and inhibit actin organization, thus impairing single cell motility, increasing collective cell migration, and reducing cancer invasion and metastasis [25].Downregulation of ADAMDEC1 was shown to upregulate active caspases 3 and 9, inhibit proliferation and induce apoptosis in glioma cells [26].We constructed a risk model using HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH expression data from a cohort of patients in TCGA, and found that it had significant prognostic value for ATRXwt glioma patients.
The tumor environment consists of cancer cells, immune cells, inflammatory cells, tumor-associated fibroblasts and various cytokines [27].Immune cells within the tumor microenvironment influence the progression of glioma [27].Glioma tissues contain a high proportion of M2 macrophages, which have the potential to promote glioma cell invasion through angiogenesis, while M1 cells have the potential to suppress angiogenesis [28,29].Natural killer cells and CD8+ T cells are susceptible to senescence in gliomas [30,31].Little has been known about the immune signature of ATRX-wt glioma, but we found that ATRX-wt patients with low risk scores tended to have higher levels of M1 and CD8+ T cells, suggesting that these infiltrating cells were able to kill cancer cells.
In order for the host to kill cancer cells, immune checkpoints need to be blocked so that deactivated cells can be reactivated.ICB therapy has shown significant curative effects in hepatocellular carcinoma and breast cancer patients [32,33].We found evidence that CD8+ T cell and M1 cell inactivation were inhibited in the tumor microenvironment of ATRX-wt glioma patients with low risk scores; thus, we analyzed whether ICB therapy was beneficial for low-risk ATRX-wt glioma patients.ICB therapy response rates were higher among low-risk patients than among high-risk patients, suggesting that ATRX-wt glioma patients with low risk scores may benefit from ICB therapy.An OncoPredict analysis predicted that high-risk glioma patients would be more sensitive to rapamycin, dasatinib, 5-fluorouracil and gemcitabine.
In conclusion, our study revealed that an immune signature based on HOXA5, PTPN2, WT1, HOXD10, POSTN, ADAMDEC1 and MYBPH expression effectively predicted the prognosis of ATRX-wt glioma patients, and demonstrated that immunotherapy was effective for low-risk patients.Our immune signature may be helpful in diagnosing and treating ATRX-wt glioma patients.

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

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.

Figure 1 .
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.

Figure 2 .
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-stromalscore 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.

Figure 3 .
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.

Figure 4 .
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.

Figure 5 .
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.

Figure 6 .
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 lowscore 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.

Figure 7 .
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 highand 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.

Figure 8 .
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.

Figure 10 .
Figure 10.Immunological characteristics of the three immune features.(A, B) Gene expression profiles of the high-and low-riskscore 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.

Figure 11 .
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.