Introduction

Gliomas are brain tumours that originate from glial cells. Star-shaped glia cells are called astrocytes, and tumours derived from these cells are referred to as astrocytomas. The 2016 WHO classification of diffuse astrocytic tumours indicates three grades with different aggressiveness1. Compared to the other types of gliomas, diffuse astrocytomas present unique challenges for treatment due to their heterogeneity, aggressive biological behaviour and diffusive growth. Diffuse astrocytomas are comprised of two major subtypes according to the newest WHO criteria: IDH-mutant and IDH-wild type (diffuse astrocytoma, anaplastic astrocytoma and glioblastoma).

The modern diagnosis of astrocytoma relies on computer-assisted tomography (CT) and magnetic resonance imaging (MRI). These techniques are indispensable for determining the localization and size of the lesions. Radiological features that can differentiate an astrocytoma from other brain parenchymal lesions have also been described2,3, but in most cases the accuracy of the radiological features is not sufficient to make clinical management decisions. Cytological and histopathological analyses of bioptates collected during brain surgery or biopsy have long been the only choice for long-term, reliable establishment of glioma grading. These considerations indicate a clear need for accurate and robust biomarkers, which may help determine the astrocytoma prognosis and grading.

An astrocytoma originates from complex interactions between developmental and genetic factors, which leads to heterogeneity across patients. To date, the transcriptomic classification divides astrocytomas into different subtypes, including classical (CL), neural (NE), mesenchymal (ME), and proneural (PN)4. Great efforts have been made to identify molecular markers for prognostic prediction. Many recent studies have focused on gene expression profiles in glioma to construct signatures5,6,7,8,9,10, which have shown great promise for prognostic prediction in individual patients. However, the limitation of the gene selection method, which uses formula designation for application within different platforms, is also obvious. Thus, identifying a more powerful and practical scoring system for the heterogeneity and intrinsic characteristics of the tumour for prognostic prediction has great clinical significance.

In the current study, with an aim of developing a scoring system for astrocytoma progression and prognosis, we carried out a supervised approach associated with the astrocytoma grade and clinical outcome. Following this principle, we performed a combined analysis to identify two robust gene sets related to astrocytoma progression and established a scoring system.

Results

Selection of astrocytoma progression-related genes and construction of the AP score

We obtained gene expression profiles for 466 astrocytomas (323 LGG and 143 GBM) and 5 NTs from TCGA. We obtained 30 genes for analysis that met the following filtering criteria (Fig.Ā 1a): (i) genes were differentially expressed between tumour and normal tissue, (ii) genes were differentially expressed between GBM and LGG tissue, (iii) genes demonstrated highly variable expression profiles (median absolute deviation, MADā€‰>ā€‰1.0), (iv) genes with the same DEG trend, and (v) genes related to patient survival (Kaplan-Meier log-rank test pā€‰<ā€‰0.05). Thus, the resulting 18 genes (POS_AP: ANXA2, CD44, DPEP1, IGBP2, IQGAP2, MMP2, MTTP, NCF1C, STEAP3, TCF19, TEAD, and TM6SF2; NEG_AP: ALOXE3, GABRD, LOC293392, LOC440905, PANX2, and SGSM1) not only represented genes that were significantly related to astrocytoma initiation and progression but also genes that exhibited heterogeneity among tumour samples and in terms of the prognosis (Fig.Ā 1a,b).

Figure 1
figure 1

Development of the AP score and its prognostic value across cohorts. (a) An overview of the AP score algorithm. The AP score algorithm uses gene expression data to output the combined degree of positive and negative regulation of astrocytoma progression. (b) Heatmap depicting the Z-score expression values of the 18 genes in the training cohort. Columns represent each sample and are labelled with their clinical characteristics, and rows represent genes and are divided into two groups representing POS_AP and NEG_AP. (c) Kaplan-Meier survival analyses based on the median cutoff AP score in the training dataset. (d) Kaplan-Meier survival analyses based on the median cutoff AP score in the GSE16011 dataset. (e) Kaplan-Meier survival analyses based on the median cutoff AP score in the REMBRANDT dataset. (f) Kaplan-Meier survival analyses based on the median cutoff AP score in TCGA GBM dataset. Mutant: IDH1 mutant, WT: IDH1 wild type; NE: Neural, PN: Pro-neural, CL: Classical, ME: Mesenchymal.

Full size image

To investigate the roles of the 18 genes in the regulation of astrocytoma progression, we nominated upregulated genes as positive regulators of astrocytoma progression (POS_AP) and downregulated genes as negative regulators of astrocytoma progression (NEG_AP).

Based on the two distinct gene sets, we employed the ssGSEA method to construct a scoring system for analysis, and the astrocytoma progression (AP) score was set as follows: AP scoreā€‰=ā€‰POS_AP score āˆ’ NEG_AP score.

The prognostic value of the AP score across different cohorts

To explore whether the AP score was related to astrocytoma patient survival, first we analysed the training dataset (TCGA RNA-seq cohort), which contained 466 cases. The samples were dichotomized into either high (nā€‰=ā€‰233) or low (nā€‰=ā€‰233) subgroups classified by the median AP score. Our analyses demonstrated that a low AP score was a potent independent marker for predicting better overall survival in the training dataset (median survival: 2835 days for the low AP score group vs 492 days for the high AP score group, Pā€‰<ā€‰0.0001, HRā€‰=ā€‰7.47; Fig.Ā 1c). We obtained similar results when analysing the GSE16011 (Pā€‰<ā€‰0.0001, HRā€‰=ā€‰2.20; Fig.Ā 1d) and REMBRANDT cohorts (Pā€‰<ā€‰0.0001, HRā€‰=ā€‰1.98; Fig.Ā 1e). These results demonstrated that the AP score was a significant reliable marker for predicting the prognosis of astrocytoma patients. To the best of our knowledge, a longitudinal dataset comprised of primary-recurrent astrocytoma is most relevant for disease progression. To validate the AP score concept, we employed the GSE4271 dataset and filtered it according to the criterion that a primary astrocytoma should be a lower-grade astrocytoma. Seven paired astrocytomas fit the criterion. Then, the AP score was calculated for the 14 astrocytomas (Fig.Ā 2a). The AP score of the recurrent group was significantly higher than that of the primary group (t test, Pā€‰=ā€‰0.02). A total of 6/7 of the tumours had an elevated AP score when recurrent (grades III to IV: 4/7; grades III to III: 2/7). The left one that showed a decreased AP score did not progress to grade IV GBM when recurrent.

Figure 2
figure 2

The AP score in the longitudinal dataset and the prognostic value of the AP score across different subgroups of the training cohort. (a) The AP score in subgroups of primary-recurrent paired astrocytomas in the GSE4271 cohort. (b) Distribution of the AP score among astrocytoma subgroups in the training cohort. ANOVA pā€‰<ā€‰0.001 in the grade and subtype subgroups. (cā€“r) Kaplan-Meier survival analyses based on the dichotomized AP score in the supergroups defined by clinicopathological factors of the training cohort.

Full size image

Distribution and prognostic value of the AP score among the astrocytoma subgroups

A related question is whether the associations between the AP score and the clinicopathological parameters and molecular features reflect the differences between glioblastoma and low-grade tumours. These conditions are so unique that any reasonable strategy will find genes that can distinguish them. To validate this question, the patients in the training cohort were stratified based on several clinicopathological factors, including age, gender, grade, received therapy, IDH mutation status, MGMT methylation status and transcriptional subtypes. We found that patients with an older age, higher grade, IDH wild type and unmethylated MGMT status demonstrated higher AP scores, whereas gender conferred little impact on the AP score distribution. With respect to transcriptional subtypes, mesenchymal and classical subtypes had relatively high AP scores, neural subtypes had the lowest AP scores, and the proneural subtype had intermediate scores (Fig.Ā 2b). Based on the median cut-off value of the AP score in the training cohort, the patients were dichotomized into either high or low supergroups within each subgroup. A high AP score accounted for a significantly high proportion of the patients with an older age, grade IV, IDH1 wild type, unmethylated MGMT status, and classical and mesenchymal subtypes, whereas a low AP score accounted for a significantly high proportion of the patients with a younger age, grades II and III, IDH1 mutant, methylated MGMT status, and neural and proneural subtypes. No significant difference was observed between the female and male patients (Supplemental TableĀ S1). We applied the dichotomized AP score for the whole cohort to all subgroups to query the prognostic value. A nearly universal result was achieved, demonstrating that a high AP score was highly correlated with a poor prognosis and vice versa (Fig.Ā 2cā€“r). With respect to the tumour grade, obtaining a significant prognostic value was difficult for grade IV cases (low AP score nā€‰=ā€‰4, high AP score nā€‰=ā€‰139). Here, a TCGA microarray dataset comprising GBM patients (all cases with wild type IDH) was used for validation (Pā€‰=ā€‰0.0363, HRā€‰=ā€‰1.27; Fig.Ā 1f). Univariate and multivariate Cox regression analyses also indicated that the AP score was an independent prognostic factor after adjusting for other clinical covariates (TableĀ 1).

Table 1 Univariate and multivariate Cox regression analyses in the training cohort.
Full size table

The AP score reliably predicted patient outcomes

Receiver operating characteristic (ROC) analysis showed that the areas under the receiver operating characteristic curve (AUCs) for the tumour grade (GBM vs lower grade astrocytoma), MGMT promoter methylation status, IDH mutation status, AP score and combined factor of the IDH status and AP score were 0.7543, 0.6509, 0.8220, 0.8367 and 0.8513, respectively, in the training cohort (Supplemental Fig.Ā S1). The results suggested that the AP score performed as well as the IDH mutation status in predicting patient outcomes and achieved an even higher prediction value when combined with the IDH mutation status.

The AP score was associated with distinct patterns of somatic mutations and copy number variation

To further investigate the impact of the AP score in the DNA level scenario, TCGA cases with available somatic mutation (460/466) and copy number variation (CNV) information (464/466) were analysed. The cases were divided into two subgroups based on the increase in the AP score (the lower and higher AP score groups).

The mutation status of well-known individual regulators of glioma was analysed (Fig.Ā 3a). Frequent mutations in IDH1, ATRX and TP53 were significantly enriched in cases with lower AP scores, whereas PTEN, EGFR, RB1, FLG, NF1, SPTA1, PIK3CA, SEMA3C, KEL, TTN, RYR2, MUC17 and PCLO were significantly enriched in cases with higher AP scores.

Figure 3
figure 3

Different mutation and copy number variation patterns of the AP score. (a) Summary of well-known individual regulators of glioma from 460 samples from the training cohort. Columns are sorted by samples with increasing AP scores. Top histogram, the sum of mutations in each sample category is indicated by the legend; Right histogram, the sum of mutations in each gene is indicated by the legend. (b) The overall copy number variation (CNV) profile in order of increasing AP score. (c,d) A distinct CNV and recurrent mutation profile is observed between astrocytomas with low and high AP scores. *Pā€‰<ā€‰0.05; **Pā€‰<ā€‰0.01; ***Pā€‰<ā€‰0.001.

Full size image

Subsequently, CNV data were investigated and showed distinct chromosomal alteration patterns between astrocytomas with low and high AP scores. Chr 7 amplification paired with Chr 10 loss was enriched in the higher AP score cases (Fig.Ā 3b). In the gene level CNV landscape, the frequently deleted genomic regions were 9p21.3 encompassing CDKN2A/CDKN2B (mean deletion, CDKN2A āˆ’0.191 vs āˆ’0.745, Pā€‰<ā€‰0.001; CDKN2B āˆ’0.191 vs āˆ’0.729, Pā€‰<ā€‰0.001) and 10q23.3 encompassing PTEN (āˆ’0.058 vs āˆ’0.635, Pā€‰<ā€‰0.001). Conversely, 7p11.2 encompassing EGFR (mean amplification, 0.195 vs 1.836, Pā€‰<ā€‰0.001), PDGFRA (4q12; 0.045 vs 0.420, Pā€‰<ā€‰0.001) and CDK4 (12q14.1; āˆ’0.029 vs 0.663, Pā€‰<ā€‰0.001) were frequently amplified with higher AP scores (Fig.Ā 3bā€“d and Supplemental Fig.Ā S2).

High AP score astrocytomas exhibited a malignant phenotype

Regarding the prognostic value of the different mutation and CNV patterns based on the AP score, GO analysis was performed to explore the functional aspects. We performed the DEG analysis based on high and low AP scores and gained 1585 genes that were upregulated in the high AP score group at a FDR of 0.01 and the lowest log fold change (logFC) of 1.5 for the GO analysis. The results suggested that a high AP score was associated with immune/inflammatory responses and cell cycle-related processes (Fig.Ā 4a). To validate the results, Pearsonā€™s correlation score (r) was calculated for each gene in the training cohort. The GO results based on 1215 genes that were positively correlated (rā€‰>ā€‰0.6) with the AP scores suggested that the genes were highly enriched in cell cycle-related processes (Fig.Ā 4b). The GO results for either positively correlated genes or differentially expressed genes were in great concordance with the cell cycle-related processes. Meanwhile, immune/inflammatory response genes were found in the DEG panel. The GO results based on the 822 downregulated genes in the high AP score group with an FDR of 0.01 and the lowest log fold change of āˆ’1.5 and the 1139 negatively correlated (rā€‰<ā€‰āˆ’0.6) with the AP score suggested that these genes were highly enriched in relatively normal nervous system functions (Supplemental Fig.Ā S3). The GSEA analysis further validated that the AP score was associated with processes or pathways that were closely related to immune/inflammatory responses and the cell cycle (Fig.Ā 3c and Supplemental Fig.Ā S3). For gene sets that were enriched in the GSEA analysis, ssGSEA was performed in the training dataset. Pearsonā€™s correlation scores were calculated based on the ssGSEA value of the gene sets queried in the GSEA and the AP score. We observed high correlations between the gene sets and AP score. The gene sets could be inferred to be of high importance in regulating astrocytoma progression. To validate this hypothesis, the samples were dichotomized into either high or low subgroups classified by the median ssGSEA score calculated from each gene set. The log-rank test results based on the classification demonstrated that all gene sets had significant prognostic value (Supplemental TableĀ S2).

Figure 4
figure 4

Biological function of the AP score. (a) The top 25 GO terms for upregulated DEGs [Log10(FC)ā€‰>ā€‰1.5] enriched based on the AP score. (b) The top 25 GO terms of highly correlated (rā€‰>ā€‰0.6) genes enriched based on the AP score. (c) GSEA results based on the increasing AP score.

Full size image

Association between the AP score and immune cells

The cell cycle and apoptosis are among the most well understood processes regulating glioma progression and the prognosis, whereas the immune/inflammatory response has only recently been identified as a key player in regulating glioma tumourigenesis and the prognosis. Thus, we performed ssGSEA with gene sets comprised of mRNA transcripts specific for most innate and adaptive immune cell subpopulations. The gene sets were representative of immune cells that constituted the ā€œimmunomeā€. The genes in each gene set were most representative of purified immune cells11. The correlation and prognostic value of each immune cell were also calculated (TableĀ 2). Macrophages, activated dendritic cells (aDCs), NK CD56dim cells, NK cells, T helper cells (Th2), eosinophils, immature dendritic cells (iDCs), neutrophils, T cells, dendritic cells (DCs), T helper cells and T helper 1 (Th1) cells, most of which represented the innate immune system, were positively correlated with the AP score and gained significant prognostic value when dichotomized by the median ssGSEA score. T central memory (Tcm) cells, T gamma delta (Tgd) cells, B cells and T effector memory (Tem) cells, most of which represented the adaptive immune system, were negatively correlated with the AP score and gained significant prognostic value. Therefore, most immune cells with significant prognostic value were associated with the AP score.

Table 2 Correlation and prognostic value of each immune cell in the training cohort.
Full size table

Discussion

As a common fatal central nervous system tumour, astrocytoma is diagnosed by genotype and histopathological criteria, although the genotype always trumps the histological phenotype1. Diffuse astrocytomas show relentless and malignant progression characterized by extensive invasion throughout the brain. Most patients with a low-grade astrocytoma progress to a highly malignant astrocytoma or glioblastoma. However, there is accumulating evidence that tumours with similar histology have distinct molecular signatures that significantly impact the treatment response and survival. Additionally, malignant advances occur within the same tumour grade, although this phenomenon has been little studied. Robust prognostic evaluation of most tumour-related factors is limited to the grade and genetic variation. Recently, the molecular classification of astrocytomas has developed rapidly. Several classification systems based on mRNA expression4,12 were established to classify tumours with equal grades into several subtypes. As tumour malignancy progresses, the tumour subtypes can transit into each other, similar to the mesenchymal and proneural subtypes in GBM13. Understanding the mechanisms of grade progression and subtype transition and blocking the major oncogenic pathways are the keys to astrocytoma therapy14.

To discover key genes regulating the oncogenesis and malignant progression of astrocytomas, we performed bioinformatics analyses of RNA expression data and revealed 18 key regulator genes corresponding to 2 subgroups, including 12 genes that were upregulated as malignancy progressed and 6 genes that were downregulated. Although the genes were mainly filtered based on GBM vs LGG and tumour vs normal tissue given other restrictions, such as same DEG trend and the prognostic analysis, we assumed that the combined signature of the two scores mostly reflected astrocytoma progression. Further functional analysis also suggested that cases with higher AP scores were associated with processes conferring tumour malignancy, whereas lower AP scores reflected relatively normal CNS function, which met the concept of astrocytoma progression.

Increasing evidence suggests that genetic changes (mutations, deletions, amplifications and overexpression) are involved in the development and progression of gliomas15. Established astrocytoma biomarkers, such as TP53 and IDH1 mutations and the recently discovered ATRX mutations, are thought to be early events in these tumours16. In our analysis, TP53, IDH1 and ATRX were more likely to be mutated in the lower AP score group than in the higher AP score group, suggesting their roles in astrocytoma oncogenesis. However, more mutation in PTEN and EGFR were detected in the high AP score group. Mutations in the tumour suppressor gene PTEN are frequent events and are associated with therapeutic resistance, because PTEN is a key player in regulating glioblastoma oncogenesis17,18. Alterations in signature oncogenes of GBM, such as EGFR, always confer a worse prognosis18,19. Different CNV patterns were also observed, such as EGFR, PDGFRA, and CDK4 amplification and CDKN2A/CDKN2B and PTEN deletion. These results suggested that the AP score could exert its progression estimation value at the DNA level.

Regarding functional aspects, a lower AP score demonstrated relatively normal CNS function, whereas a higher AP score was associated with processes conferring tumour malignancy. Among the well-studied processes regulating glioma progression and the prognosis, the cell cycle plays a central role in development and carcinogenesis20,21. Accumulated evidence indicates that abnormal cell cycle progression may confer tumour advancement and radio-resistance22 of glioma cells. Numerous studies have elucidated the important role played by the tumour microenvironment in cancer progression. In particular, the formation of hypoxic regions within the enlarged mass of solid tumours and the consequent induction of angiogenic switches are key steps in glioma development and progression23. The CNS has been clearly shown to coordinate a robust immune response with the innate and adaptive immune systems rather than inducing immune privilege24. However, the general role of the local immune response in astrocytoma progression and the prognosis remains unclear. Several studies have focused on local immune phenotypes of glioma and have indicated that either an immune signature or immune cell enrichment may play a role in predicting the patient prognosis and tumour malignancy25,26, similar to the discoveries in our study. The relatively normal CNS processes associated with a low AP score also suggest that astrocytoma formation and progression occur at the cost of sacrificing regular CNS function while gaining malignant phenotypes.

In conclusion, our findings highlighted the important role of the AP score and its related processes in the biology and clinical management of astrocytoma. The astrocytoma progression score had a considerable impact on the clinical, genomic and biological status. Evaluating the progression score of astrocytoma may help elucidate the complex role of tumour malignant processes, facilitate further astrocytoma malignancy grading and provide new insights into clinical management and drug design.

Methods

Datasets

Whole genome mRNA expression RNA-seq, somatic mutation, and copy number variation data and corresponding clinical information, including the tumour subtype, IDH mutation status, 1p-19q co-deletion status, and MGMT promoter methylation status, were downloaded from TCGA dataset16,18 (http://cancergenome.nih.gov/) as the training cohort. According to the WHO 2016 criteria, adult diffuse glioma centres around isocitrate dehydrogenase (IDH) and 1p/19q diagnostics1,27,28,29. Thus, we designed the inclusion criteria based on the newest WHO classification, with the training cohort comprised of pure 1p/19q non-codeletion cases with a mutant or wild type (WT) IDH status. The following three datasets were obtained for validation: Repository for Molecular Brain Neoplasis Data (REMBRANDT, http://caintegrator.nci.nih.gov/rembrandt), GSE16011 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?accā€‰=ā€‰GSE16011)30, and GSE4271 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?accā€‰=ā€‰gse4271)31. Additionally, we included mRNA expression microarray data from TCGA glioblastoma multiforme (GBM) project18. Normalized expression data were downloaded from each data source. Due to the lack of precise identification of the 1p/19q status in the validation cohorts, we only included those who were diagnosed as astrocytoma or glioblastoma. The training dataset comprised 466 astrocytoma and glioblastoma samples and 5 normal tissue (NT) samples. The REMBRANDT and GSE16011 datasets comprised 248 and 184 astrocytoma and glioblastoma patients, respectively, whereas the TCGA microarray dataset comprised 376 glioblastoma patients. The patientsā€™ characteristics are summarized in Supplemental TableĀ S1.

Statistical analysis

Overall survival (OS) was defined as the interval from the date of diagnosis to death or the last follow-up. The prognostic differences between patients with high or low expression of a certain gene or score (higher or lower than the median value) were calculated by the Kaplanā€“Meier method with the two-sided log-rank test using the R ā€œsurvivalā€ package. Univariate and multivariate COX regression analyses were also performed with the R ā€œsurvivalā€ package. The hazard ratio (HR) was calculated with the Mantek-Haenszel method in the GraphPad Prism 7.00 software or by the COX method in the R ā€œsurvivalā€ package using the Kaplan-Meier log-rank test or COX analysis, respectively. A two-tailed Studentā€™s t-test was performed to compare two groups of numerical values. Analysis of variance (ANOVA) was used to analyse differences among group means. The Chi-square test and Fisherā€™s exact test were used to compare frequencies between groups. The R ā€œpROCā€ package was used in the ROC analysis and for comparisons between factors. The generalized linear model was fitted when combining two factors. The prediction performance was evaluated with area under the receiver operating characteristic curve (AUC) estimation. For the ROC analysis, patients who were not censored at the last follow-up and those whose durations were less than the mean OS were excluded. Pearsonā€™s correlation analysis was used to evaluate the association between two variables. The statistical analysis was performed with the R software version 3.51 for Windows. Statistical significance was set at the level of pā€‰<ā€‰0.05.

Bioinformatics analysis

Level 3 RNA-seq data from TCGA-LGG and TCGA-GBM cases that met the inclusion criteria were downloaded and normalized within the R package ā€œTCGAbiolinksā€32. The ā€œedgeRā€ package was used to identify differentially expressed genes (DEGs) based on the threshold of a false discovery rate (FDR) less than 0.01. The ā€œTCGAbiolinksā€ package of R was applied to investigate relevant biological implications. The biological phenotype was further verified by gene set enrichment analysis (GSEA)33. Single sample GSEA (ssGSEA) of the enriched and immune cell gene sets was performed using the R gene set variation analysis (GSVA) package34. The gene list of immune cells was summarized by Gabriela Bindea et al.11.