We first explored the expression profile of members of the CXCL family. As illustrated, the mRNA expression levels of CXCL14 (P value < 0.001), CXCL9 (P value < 0.001), CXCL10 (P value < 0.001), CXCL11 (P value < 0.001), CXCL13 (P value < 0.001), CXCL3 (P value < 0.001), CXCL1 (P value < 0.001), CXCL6 (P value < 0.001) in the GBM group were significantly higher than that in the lower-grade gliomas (LGG) group, while CXCL12 (P value < 0.01) had a higher expression in the LGG group ( Figure 1A ) according to the training cohort. In the validation cohort, similar expression alternation was observed on CXCL11 (P value < 0.001), CXCL9 (P value < 0.001), CXCL10 (P value < 0.001), but not for CXCL12 and CXCL5 ( Figures 1B , C ). In LGG, the expressions of CXCL11 (P value < 0.001), CXCL9 (P value < 0.001), and CXCL10 (P value < 0.001) were higher in the Grade III group, and higher CXCL5 (P value < 0.001), CXCL2 (P value < 0.001), and CXCL3 (P value < 0.001) expression was noticed in the Grade II group ( Figure 1D ). Members like CXCL9 (P value < 0.001), CXCL10 (P value < 0.001) from the validation cohort manifested similar expression profile as that from the training cohort ( Figures 1E , F ). Members like CXCL2, CXCL3, CXCL6 in the validation cohort manifested different expression alternation comparing the results from TCGA dataset.

IDH status is a classical biomarker to predict the malignancy of gliomas, and mutant IDH gliomas showed better prognosis than IDH wild-type gliomas. Hence, we mapped the expression profile of the CXCL family based on IDH status. In the training cohort, CXCL14 (P value < 0.001), CXCL9 (P value < 0.05), CXCL10 (P value < 0.001), CXCL11 (P value < 0.001), CXCL6 (P value < 0.001), and CXCL1 (P value < 0.001) were enriched in the IDH wide-type gliomas ( Figure 1G ). Similar expression profile was also mapped in the validation cohort ( Figures 1H , I ). Together, abnormal expression of CXCLs may be associated with glioma progression.

The samples in the TCGA dataset were divided into cluster1 and cluster2 ( Supplementary Figures 1B–E ). In the TCGA dataset, samples from cluster1 showed better overall survival (OS) in gliomas (p < 0.0001) and LGG (p < 0.0001) than samples from cluster2. Nevertheless, there were no significant survival outcome differences noticed in GBM, which may constrain from its population (p = 0.12) ( Figures 2A –C ).

The support vector machine (SVM) was used to learn the characteristics of the two cluster samples and reconstruct the cluster model in validation cohort ( Supplementary Figure 1F ). Similar overall survival difference between cluster1 and cluster2 was obtained, suggesting the clustering model can predict glioma patients’ prognosis and indicated CXCL can affect glioma progression ( Figures 2D –I ).

Next, we built a risk model by employing the elastic regression analysis in order to identify the main contributor that affect glioma progression. As the correlogram showed, two co-expression clusters (CXCL1, CXCL2, CXCL3, CXLC5, and CXCL6; CXCL9, CXCL10, and CXCL11) can be mapped, implying they may share similar regulator or exert similar function ( Figure 3A ). The elastic net regression algorithm identified five members of the CXCL family as glioma prognostic-related biomarkers, and the risk of each sample was calculated according to their coefficient ( Figures 3B , C ). The expression of the CXCL family of samples of TCGA datasets, along with their clinical features, was displayed by heatmap, which was arranged by their risk score ( Figure 3D ). The expression level of CXCL9, CXCL10, CXCL11, and CXCL14 was positively correlated with risk, while CXCL12 was associated with low risk.

In the TCGA datasets, aggressiveness gliomas were more likely to be calculated with higher risk score ( Figures 3E , F ). Meanwhile, high-risk samples were also categorized as malignant glioma subtype like IDH wild-type gliomas ( Figure 3G ), MGMT unmethylated gliomas ( Figure 3H ), 1p19q non-codel gliomas ( Figure 3I ), and mesenchymal/classical gliomas ( Figure 3J ). Similar conclusion can also be obtained from the validation cohort ( Supplementary Figures 2A–J ). Therefore, the risk score of each sample associated with malignant gliomas’ clinical feature, implying its ability in predicting glioma prognosis.

Overall survival analysis suggested worse clinical outcome in high-risk samples than low-risk samples in gliomas from TCGA dataset ( Figure 4A ). In our own data, the Xiangya cohort, high-risk samples possessed shorter median survival time than low-risk samples ( Figure 4B , P = 0.0063). Similar results were also gained in other validation cohorts, including CGGA1 ( Figure 4C , P < 0.001), CGGA2 ( Figure 4D , P < 0.001), CGGA668 ( Figure 4E , P < 0.001), GSE108474 ( Figure 4F , P < 0.001), GSE43378 ( Figure 4G , P = 0.00361), GSE16011 ( Figure 4H , P < 0.001), GSE68838 ( Figure 4I , P < 0.001). In LGG, high-risk samples also have shorter median survival time than low-risk samples, whereas no significant survival outcome difference was observed in the GBM ( Supplementary Figure 3 ).

To analyze the difference between high- and low-risk samples, the GO/KEGG enrichment analyses based on the GSVA analysis were conducted. Higher enrichment score of the immune-related pathways like T cell apoptotic process, T cell–mediated immunity, T cell–mediated cytotoxicity, regulation of T cell receptor signaling pathway, antigen processing and presentation, natural killer cell–mediated cytotoxicity, and cell adhesion–associated pathways were calculated in high-risk samples ( Figures 4J ). This result implied that CXCL may be able to modulate glioma progression by affecting immunocyte function.

The expression profile of ICGs according to the risk score was mapped considering that the CXCL family plays a critical role in mediating immunocytes’ function. Positive correlation between ICG expression like HLAs, MICA, CD40LG, CD70, CD40, CTLA4, and samples’ risk was discovered ( Figure 5A ). Therefore, the microenvironment of high-risk samples may be immunosuppressed microenvironment. Then, immunocyte infiltration ratio was analyzed by conducting the ESTIMATE algorithm. The stromal score and immune score were higher in high-risk samples relative to low-risk samples, indicating the immune landscape of high-risk samples was more complicated ( Figures 5B–E , Supplementary Figures 4A, B ). Thereby, higher estimate score (the combination of stromal score and immune score) and lower tumor purity were also noticed in high-risk samples.

Immunocyte infiltration ratio was illustrated by employing the ssGSEA algorithm ( Figure 5F , Supplementary Figures 4C, D ). Results from both training cohort and validation cohort suggested that immunocytes like natural killer cells, T cells, B cells, macrophage were enriched in high-risk samples. Moreover, positive correlation between enrichment scores of immunocytes and the risk model was also noticed ( Figure 5G , Supplementary Figures 4E, F ). Together, high-risk samples infiltrated with more immunocytes and higher ICG expression, indicating its complicated, immunosuppressed microenvironment.

Temozolomide is the first-line drug for gliomas. Patients with high or low risk showed different sensitivity to temozolomide according to overall survival analysis ( Figures 6A, B , P value < 0.001). Therefore, targeting to high-risk group by combining with other compounds may improve patients’ prognosis.

Potential sensitive drugs on high-risk glioma patients are predicted as previously reported, and 17 candidate compounds were identified ( Figure 6C ). Samples of lower AUC value of compounds indicated higher sensitivity to this compound. 12-O-tetradecanoylphorbol-13-acetate, AS-703026, birinapant, CCT128930, cobimetinib, GDC-0152, LCL-161, LGX818, LY2090314, MEK162, MK2461, RITA, and Ro-4987655 are identified from the PRISM database ( Figure 6D ); Bortezomib, JW-7-52-1, and THZ-2-49 are filtered out from CTRP 1 database ( Figure 6E ). AZD5582 is filtered out from CTRP 2 database ( Figure 6F ). Notably, ro-4987655, AZD5582, and MK-2461 may serve as novel compounds to treating gliomas. Moreover, the sensitivity of compounds from PubChem database based on CellMiner website was performed ( Figures 6G , H ). NSC68516, NSC662425, NSC641205, NSC716893, NSC715229 were filtered out as high-risk glioma patients’ potential targeted compounds. Taken together, temozolomide in combination with these potential targeted compounds may slow glioma progression.

The correlation between the risk model, the cluster model, tumor grade, IDH status, 1p19q status, and MGMT status was analyzed and displayed by the Sankey diagram. The Sankey diagrams showed that glioma patients in high-risk group were related to higher-grade gliomas, IDH wild type, 1p19q non-codeletion, and cluster 2 ( Figure 7A ). Therefore, both the cluster model and the risk model can predict glioma progression. Then, the ROC curves were used to compare the prognostic ability of the risk model, the glioma histological grade, and the cluster model when taking 1p19q codel status, IDH status, and OS as different outcome variable. Results from the training cohort and validation cohort suggested that the risk model as a better prognostic predictor than the cluster model, and as efficiency as tumor pathological grade ( Figures 7B –J ).

Univariate Cox regression and multivariate Cox regression were first analyzed to identify glioma prognosis-associated clinical features. Risk (TCGA p<0.001, HR = 4.698), age (TCGA p <0.001, HR = 1.064), IDH (TCGA p <0.001, HR = 9.754), cancer (TCGA p<0.001, HR=8.750), and 1p19q (TCGA p <0.001, HR = 4.475) were considered as independent prognostic indexes for OS time ( Supplementary Figure 5 ). Therefore, those factors were further used to construct a nomogram. The global Schoenfeld test suggested risk, cancer, age, and 1p19q were qualified factors for the construction of a nomogram as previously reported ( Figures 8A –D ).

Each variable was assigned a score, and the total points can be applied to predict the patient’s survival outcome ( Figure 8E ). The accuracy of the nomogram was verified by generating the calibration curves ( Figure 8F ). The ROC curve and AUC values of this statistics in predicting 3-year and 5-year OS of gliomas patients in TCGA were 0.891 and 0.862, respectively ( Figure 8G ). In the validation cohort, it was found that the ROC and AUC of the OS predicted at 3 and 5 years were 0.814 and 0.768 in CGGA1 while 0.820 and 0.862 in CGGA2, respectively ( Figures 8H , I ). The survival outcome difference between high- and low-risk groups were observed in the training and validation cohorts ( Figures 8J –L ). Together, the nomogram constructed by the risk model showed high accuracy in predicting glioma prognosis and can be applied to clinical.

RNA‐seq data and corresponding clinical information of gliomas were obtained from the TCGA, CGGA [mRNAarray_301 (CGGA1) dataset, mRNAseq_325 (CGGA2) dataset, and mRNAseq_693 (CGGA668) dataset], and GEO database (GSE108474, GSE43378, GSE16011, and GSE68838). All expression profiles were transformed into log2(TPM+1).

For each glioma patient of the 50 samples, major exclusion criteria were incomplete follow-up data, poor quality of samples, and missing baseline clinicopathological features. Formalin-fixed paraffin-embedded tumor tissues were then collected for sequencing. One μg RNA per sample was used as input material for RNA sample preparations, and DNA was extracted and sheared followed by sequencing library preparation using NEBNext UltraTM RNA Library Prep Kit. Subsequently, PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and the Index (X) Primer. Biotin-labeled probe was applied to capture target regions after removing the PCR primer. The captured libraries were sequenced on an Illumina Hiseq platform, and 125 bp/150 bp paired-end reads were generated. Raw data (raw reads) of fastq format were first processed through in-house perlscripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, ploy-N, and low-quality reads from raw data. Meanwhile, calculation of Q20, Q30, and GC content of the clean data were performed. All downstream analyses were based on clean data with high quality. Reference genome and gene model annotation files were downloaded from the genome website directly. The reference genome index was built using Hisat2 v2.0.5, and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5, and Hisat2 was selected as the mapping tool. FeatureCounts v1.5.0-p3 was then applied to count the reads’ numbers in order to map to each gene. The TPM of each gene was calculated based on the gene length and reads count mapped to this gene. Glioma sample collection was approved by the ethics committee of Xiangya Hospital.

By using the R package “Consensus Cluster Plus” to perform consensus cluster analysis, the samples were divided into different groups to create a clustering model (44, 45). The optimum amount of Clusters was decided according to the cumulative distribution function plots and consensus matrices (46).

Support vector machine was introduced to reconstruct the cluster model in the validation cohort based on the characteristic of the cluster model with R package “e1071.” Kernel of algorithm was set as radial. Self-validation was conducted with R package “caret,” and the sensitivity of the cluster model was 0.9888 and the specificity was 0.9795 (range of sensitivity and specificity is 0 to 1).

The elastic net regression algorithm and their coefficient were calculated automatically (47). A risk system was established based on gene coefficients. According to the median of risk, patients were divided into high- and low-risk groups. Risk score is calculated as follows:

Risk = 0.028 * CXCL1 + 0.027 * CXCL9 + 0.14 * CXCL10 + 0.04 * CXCL11 + (−0.12 * CXCL12) + 0.047 * CXCL14

The GO and KEGG analyses were carried out based on GSVA analysis, and corresponding information was downloaded from the molecular signature database (MSigDB) (48, 49). Results with false discovery rate <0.05 were considered as significant (50, 51).

The ESTIMATE algorithm was applied to evaluate the composition of the tumor microenvironment (52). The immunocyte infiltration ratio was calculated by performing the ssGSEA algorithm as the previous study reported (53, 54).

Information about drug sensitivity was downloaded from the Profiling Relative Inhibition Simultaneously in Mixtures (PRISM) and the Cancer Therapeutics Response Portal (CTRP) database. Cell line expression matrix was obtained from Cancer Cell Line Encyclopedia. R package “pRRophetic” was introduced to predict sample sensitivity to certain compounds, and lower AUC values represented higher sensitivity. The knn imputation strategy was applied to imputed NA value in the expression matrix, and we discarded samples with more than 30% NA value. The “limma” package was performed to identify potential drugs, and correlation <−0.7 was set as threshold (55). Similar strategy was applied to predict drug sensitivity of compounds from Cellminer database.

Kaplan-Meier analysis was used to generate survival curves, and the validity was evaluated by log-rank test. The receiver operating characteristic (ROC) curve and the area under the curve (AUC, AUC range is 0 to 1) were introduced to compare the predictive capabilities of different models. Univariate and multivariate Cox regression analyses were used to filter prognostic variables (P‐value < 0.05). Therefore, these variables have been verified by Schoenfeld’s test to construct a nomogram with the R package “survival” and “RMS,” respectively (54, 56–59). The calibration curve and ROC were used to evaluate the accuracy of the nomogram for OS prediction.

Statistical analysis was performed by R (version 3.6.2). Wilcoxon rank‐sum test was used to compare two groups. One‐way ANOVA was used to compare multiple groups. Spearman correlation analysis was used for correlation analysis (60). NS, not statistically significant; *P < 0.05; **P < 0.01; ***P < 0.001. P‐value < 0.05 was considered as statistically significant.