Development of a T-cell activation-related module with predictive value for the prognosis and immune checkpoint blockade therapy response in glioblastoma

Patient data

A total of 672 patients with GBM and 1,151 patients with LGG were included in this study. Data on all patients were obtained from the following four databases: TCGA RNA-seq cohort (https://tcga-data.nci.nih.gov, Version 19.0), CGGA325 RNA-seq cohort (http://www.cgga.org.cn, updated on Nov 28, 2019) (Zhao et al., 2021), CGGA693 RNA-seq cohort (http://www.cgga.org.cn, updated on Nov 28, 2019) (Zhao et al., 2021), and GSE16011 microarray (updated on Mar 14, 2014). The TCGA database was used as the discovery cohort here, the CGGA325 cohort was used as the internal validation set, and CGGA693 and GSE16011 were used as the external validation sets. Samples in all datasets contained complete clinical information, such as survival, age, grade, and IDH1 mutational status. The immunotherapy dataset used was IMvigor 210 (IMvigor210CoreBiologies package), which contained 298 bladder cancer samples. The dataset included gene expression profiles, immunotherapy data, and overall survival. Overall survival (OS) was calculated from the date of diagnosis to the date of death or final follow-up. All RNA-seq expression profiles were in FPKM format, and the microarray profile was calculated by the average value of probes for each gene.

Establishment and standardization of T-cell activation-related risk score and T-cell activation priming score

The genes in the tumor-related T-cell activation gene set were subjected to univariate Cox regression analysis to obtain genes with prognostic value (Cox-P < 0.05) in the GBM cohorts (Table S1). From the intersection of TCGA GBM and CGGA325 GBM Cox results, we obtained five genes: CD276, interleukin-15 (IL-15), solute carrier family 11 member 1 (SLC11A1), TNF superfamily member 4 (TNFSF4), and triggering receptor expressed on myeloid cells like 2 (TREML2). Subsequently, the risk score was established based on the hazard ratio (HR) of these genes:

$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\Sigma [\mathrm{L}\mathrm{N}(\mathrm{H}\mathrm{R})\times \mathrm{e}\mathrm{x}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}]$$

We used the R GSVA package to get a single sample gene set enrichment analysis (ssGSEA) score of the gene set GO:0002291 (T-cell activation via T-cell receptor contact with antigen bound to MHC molecule on antigen presenting cell, http://www.broadinstitute.org/gsea/index.jsp) for each sample (Hänzelmann, Castelo & Guinney, 2013). We then obtained the T-cell activation priming score.

The min-max normalization method was used to standardize the two scores so that both would be distributed in the (0,1) interval. The specific formula was as follows:

$${\mathrm{x}}^{\ast }=(\mathrm{x}-\mathrm{m}\mathrm{i}\mathrm{n})/(\mathrm{m}\mathrm{a}\mathrm{x}-\mathrm{m}\mathrm{i}\mathrm{n})$$

x* is the standardized score, x is the original score, min is the minimum value in the dataset, and max is the maximum value in the dataset.

Establishment of immunotherapy score

We used the standardized priming score to subtract the standardized risk score. Consecutively, a Z-score standardization was performed to obtain a standard normal distribution score, which we termed immunotherapy score. The specific formula was as follows:

$$\mathrm{I}\mathrm{m}\mathrm{m}\mathrm{u}\mathrm{n}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{p}\mathrm{y}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\mathrm{Z}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(\mathrm{P}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}-\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e})$$

$$\mathrm{Z}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=(\mathrm{x}-\mu )/\delta$$

x is the original value, μ is the mean of the sample, and δ is the standard deviation of the sample.

Survival analysis

We used the R survival package and GraphPad Prism 8 to perform univariate and multivariate Cox regression analyses of patients, as well as to analyze the differences in survival between the high- and low-risk groups and the high- and low-immunotherapy score groups. We used the R survival ROC package to draw a nomogram, and then applied calibration plots to validate it.

Gene ontology analyses, kyoto encyclopedia of genes and genomes pathway analyses, and gene set enrichment analyses

We used the R limma package to calculate the differential expression of genes between high- and low-risk groups. Genes with an FDR < 0.05 were adopted and intersections were taken in the four datasets. We also used the R clusterProfiler, org.Hs.eg.db, enrichplot, and ggplot2 packages to perform GO and KEGG pathway enrichment analyses and to plot the results. We then performed GSEA analyses between the high- and low-risk groups on immune-related gene sets. Finally, we applied the R GSVA package to analyze 22 tumor-related gene-programs using the ssGSEA function.

Analyses of tumor purity and tumor immune microenvironment

The stromal score and the immune score were calculated using the R Estimate package, while the tumor purity was calculated according to the method previously described (Zhang et al., 2017). The relative quantitation of multiple immune cells in GBM was achieved in R 3.6.3 using the Metagenes method. Accordingly, immune cells were divided into three types: antitumor effect type, tumor-promoting immunosuppression type, and intermediate type (Charoentong et al., 2017). An antitumor immunity score and a protumor suppression score based on various types of immune cells were calculated for each sample. Differences in the results of the MCP counter were used to validate the above findings (Becht et al., 2016).

Tumor microenvironment immune types and hot and cold tumors

According to the level of Cytolytic activity (CYT) and PD-L1, patients were divided into four immune subtypes: TMIT I to IV (Rooney et al., 2015). In particular, TMIT I and IV subtypes were termed hot tumor, whereas TMIT II and III subtypes were termed cold tumor (Teng et al., 2015; Ock et al., 2016; Chen et al., 2017).

Statistical analysis

We mainly used R 3.6.3 and GraphPad Prism 8 for statistical analysis and plotting of the data. Univariate and multivariate Cox regression analysis were performed for survival analysis. According to the median risk score, patients were divided into high- and low-risk groups. A two-sample t-test was used to determine whether there was a significant difference between two groups. The Kaplan–Meier survival analysis and log-rank test were used to evaluate the survival difference between two groups. Calibration plots and the C-index were used to evaluate the accuracy of the nomogram model. FDR was used to evaluate the significance of GSEA results. Pearson linear correlation was used to calculate correlations, where r was the correlation coefficient. Chi-square test was used to compare the differences in immune subtypes and the sensitivity to ICB therapy between two groups (Wang et al., 2021a). In the above statistical analyses, P < 0.05 and FDR < 0.25 were regarded as statistically significant.

Establishment of T-cell activation-related risk score and clinical characteristics

As is known, the activation of T-cells in tumors has been closely related to the sensitivity of patients to immunotherapy (Kim, Herbst & Chen, 2018), and hence has been found to be heterogeneous in different samples. In order to further explore the differences in the activation of T-cells in gliomas, we obtained a total of 15 gene sets, including a T-cell activation gene set and its derivative terms from the GSEA website (http://www.gsea-msigdb.org/gsea) to be used for ssGSEA analysis (Zhu et al., 2019). Our results showed that the ssGSEA enrichment score of the T-cell activation gene set was generally higher in GBM relative to that in LGG. In addition, we found that the score was heterogeneous in the same grade gliomas (Fig. S1A).

In order to further investigate the relationship between the activation of T-cells and the clinical features in GBM, we screened out genes in GSEA T-cell activation gene set that were confirmed to be tumor-related through literature search in PubMed. And then, we obtained a tumor-related T-cell activation gene set containing 262 genes, which is summarized in Table S1. We subjected the genes to univariate Cox regression analysis, and from the intersection of TCGA GBM and CGGA325 GBM Cox results we obtained five stable prognostic genes (Cox-P < 0.05): CD276, interleukin-15 (IL-15), solute carrier family 11 member 1 (SLC11A1), TNF superfamily member 4 (TNFSF4), and triggering receptor expressed on myeloid cells like 2 (TREML2). We found that all these genes were risk indicators for poor prognosis with HR > 1, and were positively correlated to the ssGSEA enrichment score of T-cell activation-related gene sets on the GSEA database (Fig. 1A). Therefore, we constructed a prognostic risk model related to the activation of T-cells based on the sum of the LN(HR)* expression of these five genes, as follows: Risk Score = (0.5258 × CD276 expression) + (0.1477 × IL-15 expression) + (0.2178 × SLC11A1 expression) + (0.2346 × TNFSF4 expression) + (0.2192 × TREML2 expression) (Cheng et al., 2016; Zhu et al., 2019).

We observed that our generated risk score showed a significant positive correlation with the ssGSEA enrichment score (all P < 0.001, Figs. S1B–S1I), indicating that it could be a good indicator of the relative status of the activation of T-cells in GBM cohorts.

Then, we compared the risk scores on different WHO grades of gliomas in TCGA. We accordingly found that the risk score of GBMs was higher than that of LGGs (Fig. 1B). We also found that the risk score was prominently higher in the mesenchymal subtype compared with that in other subtypes in TCGA GBM (Figs. 1A, 1C). In addition, higher-risk patients were demonstrated to be mainly characterized by the presence of the IDH1 wild-type (Figs. 1A, 1D) (Venteicher et al., 2017). We subsequently confirmed all the above results in validation cohorts (Fig. S2). Our results suggest that the risk score was related to a number of clinical characteristics with poor prognosis, and a higher risk score might indicate a more complex tumor microenvironment.

Significant prognostic value of risk score

To further evaluate the prognostic value of the T-cell activation-related risk score, we divided the patients into two groups according to the median value of the risk score, and compared the differences in survival between the two groups. We used TCGA as the discovery cohort, CGGA325 as the internal validation cohort, and CGGA693 and GSE16011 as the external validation cohorts. We found that the overall survival time of patients with GBM in the high-risk group was significantly lower than that in the low-risk group. This finding was further confirmed in the validation cohorts (Fig. 2A, Figs. S3A–S3C). Moreover, we noticed that the risk score exhibited also a remarkable prognostic value in LGG and all-grade glioma (Figs. 2B–2C, Figs. S3D–S3I).

To explore whether the risk score is a continuous independent prognostic indicator, we performed univariate and multivariate Cox regression analyses. We found that the risk score was a robust indicator of poor prognosis, independent of other clinical characteristics such as grade, IDH1 mutational status, MGMT promoter methylation, radiotherapy, and age (HR = 1.605, 95% CI [1.243−2.074], P = 0.0003). We further confirmed these results in the validation sets (Figs. 2D–2E, Figs. S4A–S4F).

In order to strengthen the survival prediction and clinical translation of our risk score, we standardized the risk scores of TCGA patients with GBM using the min-max normalization method, and established a nomogram to predict the 0.5-2-year survival probability of patients with GBM (Fig. 2F). We observed that the C-index of the nomogram was 0.7455. In addition, the calibration plot showed that the predicted survival probability was highly consistent with the actual survival rate. The nomogram also showed high survival prediction accuracy in the validation sets (Fig. 2G, Figs. S4G–S4H). These results indicate that the risk score had high survival prediction accuracy and considerable clinical translation value.

High-risk score indicated a stronger immune response

To study the major differences in biological processes between the high- and low-risk groups, we obtained the differentially-expressed genes (FDR < 0.05) between the high- and low-risk groups using the limma package, and took the intersection of genes in the four datasets for enrichment analyses. The results of GO and KEGG pathway analyses reveal that immune-related biological processes (Fig. 3A, Figs. S5A–S5B, Table S2), such as inflammation, immune response regulation, interferon gamma, cell chemotaxis, lymphocyte migration, as well as some typical carcinogenic pathways, such as NF-kappa B, PI3K-AKT, P53, and apoptosis pathways were the mainly enriched processes in the high-risk group.

To further demonstrate the correlation of the risk score to certain biological processes, we performed ssGSEA on 22 canonical to cancer gene-programs (GPs). Our results showed that five gene-programs, the GP2_Immune-Tcell/Bcell, GP4_MES/ECM, GP11_Immune-IFN, GP12_Hypoxia/glycolosis, and GP15_EGF Signaling were related to the risk score in both the discovery and validation cohorts. Consistent with the above results, we found that all five gene-programs were correlated with the immune response (Fig. 3B, Figs. S5C–S5F).

Concomitantly, we performed gene set enrichment analysis (GSEA) on five immune-related terms, and found that these immune-related terms were enriched in the high-risk group (Figs. 3C–3D, Figs. S5G–S5H, Table S2). These results indicated that the high-risk group might have a more complex immune response. This finding verifies the potential differences in tumor immune microenvironments between the high- and low-risk groups.

High-risk score indicated a more complex and immune escaped tumor microenvironment

As the immune-related gene sets were demonstrated to be the mainly enriched differentially-expressed genes in the high-risk group, we investigated the relationship between the T-cell activation-related risk score and the tumor microenvironment. First, we found that the risk score showed a strongly positive correlation with tumor purity, whereas negative correlations with the immune and stromal scores. We also obtained the similar result in the validation sets (Figs. 4A–4C, Fig. S6). This finding indicated that a higher risk score was associated with lower tumor purity and a more complex tumor microenvironment.

We next explored the relationship between the risk score and various immune cellular components in the tumor microenvironment. Our Metagenes results showed that the levels of various immune cells in the tumor were increased in parallel with an increase in the risk score. According to their effects in tumors, immune cells were divided into three types: antitumor immunity type, protumor suppression type, and intermediate type. We observed that as the risk score increased, the levels of antitumor immunity and protumor suppression cells in tumors were concomitantly increased (Fig. 4D, Fig. S7) (Li et al., 2018).

We detected a significant positive correlation between the antitumor immunity and protumor suppression scores in GBM, as calculated using the Metagenes method. In particular, we found that both immune effects were enhanced with an increase in the risk score (Fig. 4E, Fig. S8). We then used the MCP counter, a tumor immune microenvironment component analysis tool, to validate our results, and found that the immune cell content was increased in tumors in parallel with an increase in the risk score (Fig. 4F, Fig. S9).

Hot tumors were mainly concentrated in the high-risk group

A previous pan cancer analysis divided tumors into four tumor microenvironment immune types based on the expression of PD-L1 and level of Cytolytic activity (CYT). Accordingly, type I and type IV with high levels of CYT were considered hot tumors, whereas type II and type III with low levels of CYT were considered cold tumors. Based on this, we initially classified GBMs according to the expression of PD-L1 and level of CYT. We found that the proportion of type I patients in the high-risk group was significantly higher than that in the low-risk group, whereas type II patients mainly in the low-risk group. We further observed that hot tumors were mainly concentrated in the high-risk group, whereas cold tumors were mainly concentrated in the low-risk group. Finally, we confirmed all these results in validation sets (Figs. 5A–5B, Fig. S10A).

Significant correlations between the risk score and expression of immune checkpoints

A number of recent studies (Tumeh et al., 2014; Du et al., 2020; Noman et al., 2020) showed that compared with cold tumors, hot tumors have increased sensitivity to immunotherapy, thus allowing for more immunotherapy opportunities. In addition, the higher expression of immune checkpoints such as PD-1 and PD-L1 has indicated that hot tumors are more suitable for targeted immunotherapy (Wang et al., 2016). We therefore calculated the correlation between the risk score and the expressions of some classic immune checkpoints, and found that the risk score was positively correlated with immune checkpoints, such as PD-1, PD-L1, CTLA-4, hepatitis A virus cellular receptor 2 (TIM3), and indoleamine 2, 3-dioxygenase 1 (IDO1) (Fig. 5C, Figs. S10B–S10D).

We inferred that the T-cell activation-related risk score might be used to evaluate the ICB response of patients with GBM (Mi et al., 2018; Anghileri et al., 2021). To this end, we used the TIDE mean (http://tide.dfci.harvard.edu/) to predict the ICB response of patients with GBM, and found that there was no difference in the proportion of patients that were sensitive to immunotherapy between the high- and low-risk groups (Fig. 5D). No FDA-approved ICB drugs for GBM are currently available. Furthermore, as no public large-scale immunotherapy clinical trial has ever been performed, we used IMvigor 210, the most widely used public immunotherapy dataset, to validate our TIDE results. Not surprisingly, we found that the T-cell activation-related risk score could not distinguish between the sensitivity of patients to immunotherapy among the two groups (Fig. 5E). Although a higher risk score was associated with an increased infiltration of immune cells in tumors, we noticed that the risk score was also significantly positively correlated with the expression of TIM-3 and IDO1 (Fig. 5C, Figs. S10B–S10D) (Wherry & Kurachi, 2015; Li et al., 2017). This finding suggests the exhaustion of most activated T-cells in tumors due to continuous stimulation by antigens (Satpathy et al., 2019). Recent studies have shown that immune checkpoint blockade therapy cannot restore the tumor-killing ability of exhausted T-cells (Yost et al., 2019). Therefore, combined with the cancer-immunity cycle (Chen & Mellman, 2013), the content of priming but not exhausted T-cells in tumors might be the key to evaluating the sensitivity of GBMs to immunotherapy.

Combining the priming and risk scores to predict the ICB response of GBMs.

The process of TCR binding to the antigen-MHC molecular complex on antigen presenting cells (APCs) is the initiation signal of the activation of T-cells (Chowell et al., 2018; Gu et al., 2021). We obtained gene signature data on the activation of T-cells via contact of their T-cell receptor with antigen bound to MHC molecules on antigen presenting cells (GO:0002291) from the GSEA database, and performed ssGSEA on GBM samples to get the priming score. Because the activation of T-cells in GBMs was significantly positively correlated with exhaustion markers, we assumed that a higher risk score in patients would be associated with greater immunosuppression in tumors (Fig. 5C, Figs. S8D–S8F). Therefore, we combined two scores to evaluate whether there were any differences in the sensitivity to ICBs among patients with GBM (Stathias et al., 2018; Wang et al., 2021b). First, we detected a positive correlation between the priming and risk scores in the discovery cohort (R² = 0.496, P < 0.001, Fig. 6A). Our TIDE prediction results revealed that the low-risk/high-priming group had the highest proportion of ICB-sensitive patients, whereas the high-risk/low-priming group had the lowest proportion of ICB-sensitive patients (P = 0.0045, chi-square test, Fig. 6B). We obtained the same results in the validation cohorts and the IMvigor 210 immunotherapy dataset (Fig. 6C, Figs. S11A–S11F). Hence, patients in the low-risk/high-priming group were demonstrated to significantly benefit from ICBs. We also observed that compared with the high-risk/low-priming patients, low-risk/high-priming patients had a longer survival in the IMvigor 210 immunotherapy cohort (Fig. 6D).

To further explore the tumor-killing potential of activated T-cells in GBMs, we standardized the priming and risk scores and then subtracted them. After further standardizing by Z-Score, we got immunotherapy score. According to the median immunotherapy score, patients were divided into a high- and low-responding group. We found that in the TCGA and validation cohorts, the TIDE predictive ICB response rate of patients with GBM with high therapy scores was significantly higher than that of patients with low therapy scores (P < 0.0001, chi-square test, Fig. 6E). Importantly, this result was also verified in validation sets (Figs. S11G–S11I). Similarly, we noticed that in the IMvigor 210 cohort, the ICB response rate in the high-immunotherapy group was significantly higher than that in the low-therapy group (P = 0.0276, chi-square test, Fig. 6F). Patients with high-therapy score had better prognostic value with longer OS (P = 0.0188, Log-rank test, Fig. 6G). Therefore, the immunotherapy score could be a good indicator of the sensitivity to ICB therapy, and patients with high immunotherapy scores would be more likely to benefit from ICBs.