The CXCL Family Contributes to Immunosuppressive Microenvironment in Gliomas and Assists in Gliomas Chemotherapy

Gliomas are a type of malignant central nervous system tumor with poor prognosis. Molecular biomarkers of gliomas can predict glioma patient’s clinical outcome, but their limitations are also emerging. C-X-C motif chemokine ligand family plays a critical role in shaping tumor immune landscape and modulating tumor progression, but its role in gliomas is elusive. In this work, samples of TCGA were treated as the training cohort, and as for validation cohort, two CGGA datasets, four datasets from GEO database, and our own clinical samples were enrolled. Consensus clustering analysis was first introduced to classify samples based on CXCL expression profile, and the support vector machine was applied to construct the cluster model in validation cohort based on training cohort. Next, the elastic net analysis was applied to calculate the risk score of each sample based on CXCL expression. High-risk samples associated with more malignant clinical features, worse survival outcome, and more complicated immune landscape than low-risk samples. Besides, higher immune checkpoint gene expression was also noticed in high-risk samples, suggesting CXCL may participate in tumor evasion from immune surveillance. Notably, high-risk samples also manifested higher chemotherapy resistance than low-risk samples. Therefore, we predicted potential compounds that target high-risk samples. Two novel drugs, LCL-161 and ADZ5582, were firstly identified as gliomas’ potential compounds, and five compounds from PubChem database were filtered out. Taken together, we constructed a prognostic model based on CXCL expression, and predicted that CXCL may affect tumor progression by modulating tumor immune landscape and tumor immune escape. Novel potential compounds were also proposed, which may improve malignant glioma prognosis.

Introduction

A glioma, accounting for 81% of malignant brain tumors, is a primary brain tumor originating from glial stem cells or progenitor cells (1, 2). The malignancy of gliomas is attributed to its rapid cell proliferation and abnormal angiogenesis (3). World Health Organization graded gliomas from I to IV based on tumor histological features. Grade IV gliomas, also known as GBM, progress aggressively, show high recurrence rate, and are resistant to tumor treatment, along with median survival time less than 14.6 months. Molecular biomarkers like IDH status, MGMT promoter status, 1p19q, ARTX have been considered as glioma progression-associated marker. Molecular classification like mesenchymal, classical, and proneural was also proposed to assist in predicting glioma progression. Current standard treatment for gliomas includes maximal surgical removal combined with radio-chemotherapy, but patients’ survival outcome is still unsatisfactory (4). Therefore, the exploration of potential biomarkers may benefit glioma patients and assist in clinical treatment decision.

Tumor microenvironment (TME) consists of tumor cells and non-tumor cells, such as microglia, peripheral macrophages, myeloid-derived suppressor cells, vascular endothelial cells, and others (5, 6). The infiltration of those non-tumor cells has been proven to modulate tumor progression and tumor treatment sensitivity. Therefore, chemokines secreted by gliomas affect not only TME components but also tumor progression (7).

The C-X-C motif chemokine ligand (CXCL) family widely participates in immunocyte recruitment and affecting tumor progression like tumor migration and angiogenesis (8, 9). For instance, CXCL1 promotes tumor angiogenesis in ovarian cancer (10). CXCL1, CXCL2, and CXCL8 participate in the formation of endothelial tube in cervical cancer (11). CXCL1, CXCL5, and CXCL16 are able to facilitate tumor metastasis like lung cancer (12) and gastric cancer (13). In gliomas, abnormal expression profiles of CXCL9, CXCL10, CXCL11, and CXCL12 were noticed among different pathological grade gliomas, implying their potential relationship with glioma progression (14–16). Nevertheless, their specific role in gliomas is elusive.

In our study, samples of TCGA were set as training cohort, while CGGA2, CGGA1, four datasets from GEO database, and our own clinical samples were treated as validation cohort (Supplementary Figure 1A). The expression profile of the CXCL family was first depicted. Then, consensus clustering analysis was introduced to establish the cluster model based on CXCL expression profile. Abnormal CXCL expression was noticed within the cluster model implying their role in glioma progression. In order to identify the main contributor of glioma progression and enhance the cluster model prognosis prediction accuracy, the elastic net analysis was introduced and calculated risk score of each sample. CXCL9, CXCL10, CXCL11, CXCL12, and CXCL14 were filtered as main contributors of glioma progression. High-risk samples exhibited worse clinical outcome, indicating the prognostic prediction ability of the risk model. Besides, high immunocyte infiltration ratio and immune checkpoint gene (ICG) expression were also discovered in high-risk samples. Therefore, potential compounds targeted to high-risk samples were also predicted.

Results

Abnormal Expression Profile of the CXCL Family May Affect Glioma Progression

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.

FIGURE 1

Figure 1 Relationships between CXCL expression patterns and clinical characters of gliomas. The expression profile of the CXCL family in the LGGGBM group (A) from the TCGA database and the validation cohort (B, C). The heatmap of the CXCL family expression in LGG group (D–F) and GBM group (G–I) from the training and validation cohorts. *p < 0.05, **p < 0.01, ***p < 0.001.

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.

Cluster2 Samples Exhibit Aggressiveness Growth Pattern Than Cluster1 Samples

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

FIGURE 2

Figure 2 The prognostic value of the cluster model in gliomas. Kaplan–Meier survival analysis were used to show overall survival outcome difference between the two clusters in LGG, GBM, and LGGGBM samples from TCGA (A–C), CGGA1 (D–F), and CGGA2 (G–I).

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

Constructing the Risk Model Based on Elastic Net Regression Analysis

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.

FIGURE 3

Figure 3 Constructing the risk model. (A) The co-expression network of the CXCL family. (B, C) The construction of the risk model based on the expression profile of the CXCL family by performing elastic net regression algorithm. (D) Heatmap displayed the alternation of the CXCL family members’ expression according to the risk model, and corresponding clinical features were also mapped. The distribution of the risk in gliomas’ pathological grade (E), IDH status (F), cancer type (G), MGMG status (H), 1p19q status (I), and subtype (J). NS, not significant, ***p < 0.001.

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

FIGURE 4

Figure 4 The overall survival analysis and biofunction prediction based on the risk model. Kaplan–Meier survival analysis based on the risk model from the training cohort (A) and the validation cohort, including Xiangya cohort (B), CGGA1 (C), CGGA2 (D), CGGA688 (E), GSE108474 (F), GSE43378 (G), GSE16011 (H), GSE68838 (I). (J) GO/KEGG enrichment analysis based on the GSVA analysis in the training cohort.

High-Risk Samples Showed Higher Immunocyte Infiltration Ratio and ICGs Expression

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.

FIGURE 5

Figure 5 The expression of ICGs and immunocyte infiltration ratio. (A) ICG expression was mapped based on the risk model. The ImmuneScore (B), StromalScore (C), ESTIMATEscore (D), and purity (E) difference between high- and low-risk group in TCGA database. (F) Immunocyte infiltration ratio enrichment score based on the risk model. (G) The correlation between the enrichment score of immunocytes and the risk model. NS, not significant, ***p < 0.001.

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.

Chemotherapy Suggestion Based on the Risk Model

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.

FIGURE 6

Figure 6 Chemotherapy suggestion based on the risk model. Chemotherapy efficacy difference between high- and low-risk group in the CGGA1 (A) and CGGA2 (B) datasets. (C) Correlation between risk and the AUC value of 17 candidate compounds. (D–F) Distribution of the AUC value of candidate compounds in the risk model. (G) Correlation between compounds’ IC50 and risk from Cellminer dataset. (H) The difference between IC50 of compounds and the risk model ***p < 0.001.

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 Application of the Risk Model in Clinical

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

FIGURE 7

Figure 7 Comparing the prognostic ability of the cluster model, the risk model, and glioma pathological grade. (A) The Sankey diagram revealed the potential connection between glioma pathological grade, risk, cluster, IDH status, 1p19q status, and MGMT status. ROC curve generated based on the risk model by taking the IDH status (B), OS (C), and 1p19q status (D) as outcome variable in the training cohort. The validation of ROC curve in the CGGA1 (E–G) and CGGA2 (H–J) database.

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

FIGURE 8

Figure 8 Prognostic nomogram based on the risk model. (A–D) The Schoenfeld test of the factors involved in the construction of the nomogram. (E) The nomogram based on the risk model. (F) The calibration curve of 3-year and 5-year OS based on the nomogram. ROC curves and AUC values from the nomogram of 3-year and 5-year OS, in TCGA datasets (G), CGGA1 datasets (H), and CGGA2 datasets (I). (J–L) Survival analysis based on the nomogram in TCGA datasets (J), CGGA1 datasets (K), and CGGA2 datasets (L).

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.

Discussion

Multiple studies reported that the CXCL family played a critical role in tumorigenesis, tumor cell proliferation and metastasis, tumor cell resistance to drugs, and tumor angiogenesis (17–20). For instance, CXCL14 promoted GBM progression by modulating tumor cell proliferation and migration (21). The STAT3 inhibitor can target the central nervous system tumor by induing immunocyte tumor homing in a CXCL10-dependent manner (22). CXCL12 was valued as a prognostic biomarker of LGG (23). In this work, CXCL9, CXCL10, and CXCL14 were also identified as potential vital regulators of glioma progression among the CXCL family.

The immune evasion mechanism of gliomas plays an important role in glioma tumor resistance to treatment and tumorigenesis (24–26). Biofunction prediction suggested that the CXCL family may promote glioma progression through inducing immune escape and affecting immunocyte infiltration. CD70 is a critical mediator of immunocytes’ activation in the tumor microenvironment (27, 28). High expression of classical ICGs like CTLA4 was also noticed in high-risk glioma samples, implying its immunosuppressive microenvironment (29).

Immunocytes like CD4 T cell, CD8 T cell, activated T cell, and memory T cell were preferentially infiltrated in high-risk samples. High immunocyte infiltration usually suggests high immunogenicity. However, a previous study reported that high ICG expression and function-impaired T cells in GBM together contribute to an immunosuppressed tumor microenvironment (30). Therefore, ICGs interfered in immunocytes’ function in high-risk samples, in the end resulting in tumor evasion from immune surveillance. For instance, PD-1 can inhibit immunocytes’ biofunction in gliomas (29, 31). CTLA-4 facilitates immunosuppressive microenvironment by inhibiting antigen-specific T cell activation and enhancing myeloid-derived suppressor cells (32, 33). Together, considering the CXCL family can induce immunosuppressive microenvironment and immunocyte infiltration in solid tumor (34–36), targeting the members of CXCL family may improve tumor immunotherapy efficacy.

A previous study reported that the CXCL family participated in gliomas’ resistance to chemotherapy (37), and similar conclusion can be obtained from our study. Therefore, novel compounds targeted to high-risk samples may improve that situation. Novel high-risk sample-targeted therapeutic drugs like Ro-4987655, AZD5582, and MK-2461 have been reported to play a role in other tumors, but not in gliomas. For example, Ro-4987655, a highly selective mitogen-activated protein kinase kinase (MEK) inhibitor (38), has proved its efficacy in BRAF V600 mutant melanoma, BRAF wild-type melanoma, and KRAS mutant NSCLC patients (39). AZD5582 can trigger cell apoptosis in pancreatic cancer cells (40) and non-small-cell lung cancer (41). MK-2461 is an ATP-competitive multitargeted inhibitor of activated c-Met and able to slow the progression of pancreatic cancer (42, 43). Additionally, compounds from PubChem database, like NSC168516 and NSC715229, were also predicted as high-risk glioma sample–sensitive drugs. Taken together, glioma patients’ survival time may be able to be prolonged by combining those compounds with temozolomide.

In conclusion, our research revealed the expression characteristics of the CXCL gene family and constructed a high-precision glioma prognosis model. Moreover, this model highlighted the relationship between CXCL and tumor immunogenicity and offered novel treatment strategies.

Materials and Methods

Data Processing

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.

Construction of Prognostic Model

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

Biological Function Prediction

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

Potential Compounds Prediction

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.

Survival Analysis and Nomogram

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

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.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author Contributions

ZW, YL, and YM: Writing—original draft, writing—review, editing, data curation, and formal analysis. HZ and ZD: Investigation and methodology. WY, HC, and XZ: Validation. ZL and QC: Conceptualization, supervision, project administration, and funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

This study is supported by the National Nature Science Foundation of China (NO. 82073893, NO. 81873635, NO. 81703622); the China Postdoctoral Science Foundation (NO.2018M633002); the Natural Science Foundation of Hunan Province (NO. 2018JJ3838, NO. 2018SK2101); the Hunan Provincial Health and Health Committee Foundation of China (C2019186); and Xiangya Hospital Central South University postdoctoral foundation, and the Fundamental Research Funds for the Central Universities of Central South University (No. 2021zzts1027).

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

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

References

1. Weller M, Wick W, Aldape K, Brada M, Berger M, Pfister SM, et al. Glioma. Nat Rev Dis Primers (2015) 1:15017. doi: 10.1038/nrdp.2015.17

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ostrom QT, Bauchet L, Davis FG, Deltour I, Fisher JL, Langer CE, et al. The Epidemiology of Glioma in Adults: A "State of the Science" Review. Neuro Oncol (2014) 16:896–913. doi: 10.1093/neuonc/nou087

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Peng Z, Liu C, Wu M. New Insights Into Long Noncoding RNAs and Their Roles in Glioma. Mol Cancer (2018) 17:61. doi: 10.1186/s12943-018-0812-2

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Smits M, van den Bent MJ. Imaging Correlates of Adult Glioma Genotypes. Radiology (2017) 284:316–31. doi: 10.1148/radiol.2017151930

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Yan X, Xie Y, Yang F, Hua Y, Zeng T, Sun C, et al. Comprehensive Description of the Current Breast Cancer Microenvironment Advancements via Single-Cell Analysis. J ExpClin Cancer Res (2021) 40(1):142. doi: 10.1186/s13046-021-01949-z

CrossRef Full Text | Google Scholar

6. Gieryng A, Pszczolkowska D, Walentynowicz KA, Rajan WD, Kaminska B. Immune Microenvironment of Gliomas. Lab Invest (2017) 97:498–518. doi: 10.1038/labinvest.2017.19

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Chen Z, Hambardzumyan D. Immune Microenvironment in Glioblastoma Subtypes. Front Immunol (2018) 9:1004. doi: 10.3389/fimmu.2018.01004

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Vandercappellen J, Van Damme J, Struyf S. The Role of CXC Chemokines and Their Receptors in Cancer. Cancer Lett (2008) 267:226–44. doi: 10.1016/j.canlet.2008.04.050

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Wu L, Awaji M, Saxena S, Varney ML, Sharma B, Singh RK. IL-17-CXC Chemokine Receptor 2 Axis Facilitates Breast Cancer Progression by Up-Regulating Neutrophil Recruitment. Am J Pathol (2020) 190:222–33. doi: 10.1016/j.ajpath.2019.09.016

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ouh YT, Cho HW, Lee JK, Choi SH, Choi HJ, Hong JH. CXC Chemokine Ligand 1 Mediates Adiponectin-Induced Angiogenesis in Ovarian Cancer. Tumour Biol (2019) 42:1010428319842699. doi: 10.1177/1010428319842699

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zhang W, Wu Q, Wang C, Yang L, Liu P, Ma C. AKIP1 Promotes Angiogenesis and Tumor Growth by Upregulating CXC-Chemokines in Cervical Cancer Cells. Mol Cell Biochem (2018) 448:311–20. doi: 10.1007/s11010-018-3335-7

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Liang K, Liu Y, Eer D, Liu J, Yang F, Hu K. High CXC Chemokine Ligand 16 (CXCL16) Expression Promotes Proliferation and Metastasis of Lung Cancer via Regulating the NF-kappaB Pathway. Med Sci Monit (2018) 24:405–11. doi: 10.12659/MSM.906230

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zhou Z, Xia G, Xiang Z, Liu M, Wei Z, Yan J, et al. A C-X-C Chemokine Receptor Type 2-Dominated Cross-Talk Between Tumor Cells and Macrophages Drives Gastric Cancer Metastasis. Clin Cancer Res (2019) 25:3317–28. doi: 10.1158/1078-0432.CCR-18-3567

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Weenink B, Draaisma K, Ooi HZ, Kros JM, Sillevis Smitt PAE, 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

PubMed Abstract | CrossRef Full Text | Google Scholar

15. do Carmo A, Patricio I, Cruz MT, Carvalheiro H, Oliveira CR, Lopes MC. CXCL12/CXCR4 Promotes Motility and Proliferation of Glioma Cells. Cancer Biol Ther (2010) 9:56–65. doi: 10.4161/cbt.9.1.10342

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Rao S, Sengupta R, Choe EJ, Woerner BM, Jackson E, Sun T, et al. CXCL12 Mediates Trophic Interactions Between Endothelial and Tumor Cells in Glioblastoma. PloS One (2012) 7:e33005. doi: 10.1371/journal.pone.0033005

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Lukaszewicz-Zajac M, Paczek S, Muszynski P, Kozlowski M, Mroczko B. Comparison Between Clinical Significance of Serum CXCL-8 and Classical Tumor Markers in Oesophageal Cancer (OC) Patients. Clin Exp Med (2019) 19:191–9. doi: 10.1007/s10238-019-00548-9

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Zhang G, Luo X, Zhang W, Chen E, Xu J, Wang F, et al. CXCL-13 Regulates Resistance to 5-Fluorouracil in Colorectal Cancer. Cancer Res Treat (2020) 52:622–33. doi: 10.4143/crt.2019.593

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Rodrigues G, Hoshino A, Kenific CM, Matei IR, Steiner L, Freitas D, et al. Tumour Exosomal CEMIP Protein Promotes Cancer Cell Colonization in Brain Metastasis. Nat Cell Biol (2019) 21:1403–12. doi: 10.1038/s41556-019-0404-4

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zhang Z, Tan X, Luo J, Cui B, Lei S, Si Z, et al. GNA13 Promotes Tumor Growth and Angiogenesis by Upregulating CXC Chemokines via the NF-kappaB Signaling Pathway in Colorectal Cancer Cells. Cancer Med (2018) 7:5611–20. doi: 10.1002/cam4.1783

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Fazi B, Proserpio C, Galardi S, Annesi F, Cola M, Mangiola A, et al. The Expression of the Chemokine CXCL14 Correlates With Several Aggressive Aspects of Glioblastoma and Promotes Key Properties of Glioblastoma Cells. Int J Mol Sci (2019) 20. doi: 10.3390/ijms20102496

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Fujita M, Zhu X, Sasaki K, Ueda R, Low KL, Pollack IF, et al. Inhibition of STAT3 Promotes the Efficacy of Adoptive Transfer Therapy Using Type-1 CTLs by Modulation of the Immunological Microenvironment in a Murine Intracranial Glioma. J Immunol (2008) 180:2089–98. doi: 10.4049/jimmunol.180.4.2089

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Salmaggi A, Gelati M, Pollo B, Marras C, Silvani A, Balestrini MR, et al. CXCL12 Expression Is Predictive of a Shorter Time to Tumor Progression in Low-Grade Glioma: A Single-Institution Study in 50 Patients. J Neurooncol (2005) 74:287–93. doi: 10.1007/s11060-004-7327-y

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Fehervari Z. Glioma Immune Evasion. Nat Immunol (2017) 18:487. doi: 10.1038/ni.3736

CrossRef Full Text | Google Scholar

25. Daubon T, Hemadou A, Romero Garmendia I, Saleh M. Glioblastoma Immune Landscape and the Potential of New Immunotherapies. Front Immunol (2020) 11:585616. doi: 10.3389/fimmu.2020.585616

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kane JR, Zhao J, Tsujiuchi T, Laffleur B, Arrieta VA, Mahajan A, et al. CD8(+) T-Cell-Mediated Immunoediting Influences Genomic Evolution and Immune Evasion in Murine Gliomas. Clin Cancer Res (2020) 26:4390–401. doi: 10.1158/1078-0432.CCR-19-3104

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Jin L, Ge H, Long Y, Yang C, Chang YE, Mu L, et al. CD70, a Novel Target of CAR T-Cell Therapy for Gliomas. Neuro Oncol (2018) 20:55–65. doi: 10.1093/neuonc/nox116

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Zhao J, Que W, Du X, Fujino M, Ichimaru N, Ueta H, et al. Monotherapy With Anti-CD70 Antibody Causes Long-Term Mouse Cardiac Allograft Acceptance With Induction of Tolerogenic Dendritic Cells. Front Immunol (2020) 11:555996. doi: 10.3389/fimmu.2020.555996

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang H, Dai Z, Wu W, Wang Z, Zhang N, Zhang L, et al. Regulatory Mechanisms of Immune Checkpoints PD-L1 and CTLA-4 in Cancer. J Exp Clin Cancer Res (2021) 40:184. doi: 10.1186/s13046-021-01987-7

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Woroniecka K, Chongsathidkiet P, Rhodin K, Kemeny H, Dechant C, Farber SH, et al. T-Cell Exhaustion Signatures Vary With Tumor Type and Are Severe in Glioblastoma. Clin Cancer Res (2018) 24:4175–86. doi: 10.1158/1078-0432.CCR-17-1846

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Galstyan A, Markman JL, Shatalova ES, Chiechi A, Korman AJ, Patil R, et al. Blood-Brain Barrier Permeable Nano Immunoconjugates Induce Local Immune Responses for Glioma Therapy. Nat Commun (2019) 10:3850. doi: 10.1038/s41467-019-11719-3

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Huang J, Liu F, Liu Z, Tang H, Wu H, Gong Q, et al. Immune Checkpoint in Glioblastoma: Promising and Challenging. Front Pharmacol (2017) 8:242. doi: 10.3389/fphar.2017.00242

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Liu F, Huang J, Liu X, Cheng Q, Luo C, Liu Z. CTLA-4 Correlates With Immune and Clinical Characteristics of Glioma. Cancer Cell Int (2020) 20:7. doi: 10.1186/s12935-019-1085-6

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Che Y, Yang Y, Suo J, An Y, Wang X. Induction of Systemic Immune Responses and Reversion of Immunosuppression in the Tumor Microenvironment by a Therapeutic Vaccine for Cervical Cancer. Cancer Immunol Immunother (2020) 69:2651–64. doi: 10.1007/s00262-020-02651-3

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Gonzalez-Aparicio M, Alfaro C. Significance of the IL-8 Pathway for Immunotherapy. Hum Vaccin Immunother (2020) 16:2312–7. doi: 10.1080/21645515.2019.1696075

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Farsaci B, Donahue RN, Coplin MA, Grenga I, Lepone LM, Molinolo AA, et al. Immune Consequences of Decreasing Tumor Vasculature With Antiangiogenic Tyrosine Kinase Inhibitors in Combination With Therapeutic Vaccines. Cancer Immunol Res (2014) 2:1090–102. doi: 10.1158/2326-6066.CIR-14-0076

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Bruyere C, Mijatovic T, Lonez C, Spiegl-Kreinecker S, Berger W, Kast RE, et al. Temozolomide-Induced Modification of the CXC Chemokine Network in Experimental Gliomas. Int J Oncol (2011) 38:1453–64.

PubMed Abstract | Google Scholar

38. Akinleye A, Furqan M, Mukhi N, Ravella P, Liu D. MEK and the Inhibitors: From Bench to Bedside. J Hematol Oncol (2013) 6:27. doi: 10.1186/1756-8722-6-27

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zimmer L, Barlesi F, Martinez-Garcia M, Dieras V, Schellens JH, Spano JP, et al. Phase I Expansion and Pharmacodynamic Study of the Oral MEK Inhibitor RO4987655 (CH4987655) in Selected Patients With Advanced Cancer With RAS-RAF Mutations. Clin Cancer Res (2014) 20:4251–61. doi: 10.1158/1078-0432.CCR-14-0341

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Moon JH, Shin JS, Hong SW, Jung SA, Hwang IY, Kim JH, et al. A Novel Small-Molecule IAP Antagonist, AZD5582, Draws Mcl-1 Down-Regulation for Induction of Apoptosis Through Targeting of Ciap1 and XIAP in Human Pancreatic Cancer. Oncotarget (2015) 6:26895–908. doi: 10.18632/oncotarget.4822

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Hao Q, Tang H. Interferon-Gamma and Smac Mimetics Synergize to Induce Apoptosis of Lung Cancer Cells in a TNFalpha-Independent Manner. Cancer Cell Int (2018) 18:84. doi: 10.1186/s12935-018-0579-y

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Pan BS, Chan GK, Chenard M, Chi A, Davis LJ, Deshmukh SV, et al. MK-2461, a Novel Multitargeted Kinase Inhibitor, Preferentially Inhibits the Activated C-Met Receptor. Cancer Res (2010) 70:1524–33. doi: 10.1158/0008-5472.CAN-09-2541

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Inoue K, Ohtsuka H, Tachikawa M, Motoi F, Shijo M, Douchi D, et al. MK2461, a Multitargeted Kinase Inhibitor, Suppresses the Progression of Pancreatic Cancer by Disrupting the Interaction Between Pancreatic Cancer Cells and Stellate Cells. Pancreas (2017) 46:557–66. doi: 10.1097/MPA.0000000000000778

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wilkerson MD, Hayes DN. ConsensusClusterPlus: A Class Discovery Tool With Confidence Assessments and Item Tracking. Bioinformatics (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Datta S, Datta S. Comparisons and Validation of Statistical Clustering Techniques for Microarray Gene Expression Data. Bioinformatics (2003) 19:459–66. doi: 10.1093/bioinformatics/btg025

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Goeman JJ. L1 Penalized Estimation in the Cox Proportional Hazards Model. Biom J (2010) 52:70–84.

PubMed Abstract | Google Scholar

48. Hanzelmann S, Castelo R, Guinney J. GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

49. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) Hallmark Gene Set Collection. Cell Syst (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Frost HR, Shen L, Saykin AJ, Williams SM, Moore JH, Alzheimer's Disease Neuroimaging I. Identifying Significant Gene-Environment Interactions Using a Combination of Screening Testing and Hierarchical False Discovery Rate Control. Genet Epidemiol (2016) 40:544–57. doi: 10.1002/gepi.21997

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Miller RE, Breheny P. Marginal False Discovery Rate Control for Likelihood-Based Penalized Regression Models. Biom J (2019) 61:889–901. doi: 10.1002/bimj.201800138

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-Cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Wang Z, Su G, Dai Z, Meng M, Zhang H, Fan F, et al. Circadian Clock Genes Promote Glioma Progression by Affecting Tumour Immune Infiltration and Tumour Cell Proliferation. Cell Prolif (2021) 54:e12988. doi: 10.1111/cpr.12988

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Yang C, Huang X, Li Y, Chen J, Lv Y, Dai S. Prognosis and Personalized Treatment Prediction in TP53-Mutant Hepatocellular Carcinoma: An In Silico Strategy Towards Precision Oncology. Brief Bioinform (2021) 22.

Google Scholar

56. Nunez E, Steyerberg EW, Nunez J. Regression Modeling Strategies. Rev Esp Cardiol (2011) 64:501–7.

PubMed Abstract | Google Scholar

57. Nick TG, Hardin JM. Regression Modeling Strategies: An Illustrative Case Study From Medical Rehabilitation Outcomes Research. Am J Occup Ther (1999) 53:459–70. doi: 10.5014/ajot.53.5.459

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Schulgen G, Olschewski M, Krane V, Wanner C, Ruf G, Schumacher M. Sample Sizes for Clinical Trials With Time-to-Event Endpoints and Competing Risks. Contemp Clin Trials (2005) 26:386–96. doi: 10.1016/j.cct.2005.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Zhang N, Dai Z, Wu W, Wang Z, Cao H, Zhang Y, et al. The Predictive Value of Monocytes in Immune Microenvironment and Prognosis of Glioma Patients Based on Machine Learning. Front Immunol (2021) 12:656541. doi: 10.3389/fimmu.2021.656541

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Qiu H, Tian W, He Y, Li J, He C, Li Y, et al. Integrated Analysis Reveals Prognostic Value and Immune Correlates of CD86 Expression in Lower Grade Glioma. Front Oncol (2021) 11:654350. doi: 10.3389/fonc.2021.654350

PubMed Abstract | CrossRef Full Text | Google Scholar