Prognostic impact of a lymphocyte activation-associated gene signature in GBM based on transcriptome analysis

Datasets description

We downloaded three gene expression datasets (TCGA, GSE16011 (Gravendeel et al. , 2009) and GSE13041 (Lee, Scheck & Nelson, 2008)) from The Cancer Genome Atlas (TCGA) portal (https://portal.gdc.cancer.gov) and Gene Expression Omnibus (GEO) (Edgar, Domrachev & Lash, 2002; Barrett et al., 2013), which contained patient outcome and clinicopathological factors. The gene expression profile of each dataset was normalized by replacing the expression level e with log2(e+1). The patients without matched clinical information were excluded. Consequently, a total of 866 GBM patients (525 GBM patients from TCGA, 191 GBM patients from GSE13041 and 150 GBM patients from GSE16011) were included in this study. The TCGA dataset was treated as a training dataset, and another two datasets (GSE16011 and GSE13041) were used as independent validation cohorts. In the training dataset, the median survival time of GBM patients was 12.4 months (ranging from 0.1 to 127.5 months). And there were 144 classical, 155 mesenchymal, 83 neural, and 99 proneural GBM samples. In Table 1 we summarized clinicopathological data belonging to GBM patients in the TCGA dataset. The median survival of GBM patients was 13.0 months (ranging from 0.2 to 111.8 months) and 8.7 months (ranging from 0.2 to 150.7 months) in the validation datasets, respectively (Table 2).

Identifying the lymphocyte activation-associated gene signature

Differential expression analysis was performed based on transcriptome profile of GBM patients by the “limma” R package in the training dataset. Genes with the cutoff criteria of —log2-fold change—≥1 and FDR < 0.000001 between tumor and normal tissues were regarded as differentially expressed genes (DEGs). Using the univariate Cox proportional hazards regression algorithm, we retained the DEGs that were strongly associated with the survival of patients with GBM (P < 0.05). We focused on the prognostic genes that participated in the lymphocyte activation process. The lymphocyte activation-related genes were collected by GO terms in Gene Ontology (Ashburner et al., 2000; Gene Ontology, 2021). Then, we further selected genes associated with overall survival by stepwise regression analysis in the training dataset. Ultimately, we kept genes with the smallest Akaike information criteria (AIC) value to construct the prognostic model. There were four genes (TCF3, IGFBP2, TYRO3 and NOD2) that were identified in our study (Fig. S1). Finally, a risk score model was developed based on the gene expression weighted by regression coefficients of univariable Cox regression: Risk score = (−0.6246210 × expression level of TCF3) + (0.2992  × expression level of IGFBP2) + (0.2421068  × expression level of TYRO3) + (0.2155469  × expression level of NOD2). Based on this model, the patients were classified into the high- and low-risk groups by the median risk score in the training dataset.

Statistical analysis

In the survival analysis, we analyzed only overall survival as the end point. The survival differences between the two groups were visualized by Kaplan–Meier analysis and were compared by the log-rank test. The prognostic value of the lymphocyte activation-associated gene signature was estimated by univariate and multivariate Cox proportional hazards regression models. The Cox proportional hazards regression model was used to calculate hazard ratio (HR) and 95% confidence intervals (CI). The P-values smaller than 0.05 were considered to be statistically significant. We obtained gene signatures (that is cell markers) of lymphocytes (including natural killer cells, T cells and B cells) from the CIBERSORT web portal (Newman et al., 2015) and MCP-counter web portal (Becht et al., 2016). We combined information from these two sources to get gene signatures of lymphocytes. The lymphocyte activity of each patient was calculated by single sample gene set enrichment analysis (ssGSEA) (Barbie et al., 2009) based on the gene signatures of lymphocytes. Each ssGSEA enrichment score represents the degree to which the genes in a particular gene set are coordinately regulated within a sample. In our study, the gene signature of lymphocytes means the gene set of ssGSEA and the lymphocyte activity means the enrichment score. The concordance index (C-index) was used to compare the prognostic efficacy among this signature, age, sex, and the combined model using a logistic regression with the aforementioned three variables, and the significant p-value were calculated by rcorrp.cens function in Hmisc package. R software (http://www.r-project.org) was used to perform all of the statistical analyses, version 3.5.1 (Packages: survival, survminer, limma, Hmisc, ggplot2, GSVA).

Development of a prognostic signature related to lymphocyte activation

The gene expression profile and clinical information of GBM patients were obtained from TCGA, which were treated as the training dataset. There were 525 tumor samples and 10 normal samples. The tumor and normal samples came from different samples. We identified 2,364 differentially expressed genes between tumor and normal samples by using the limma package (FDR < 0.000001 and —log2-fold change— ≥1), including 933 up-regulated and 1,431 down-regulated genes (Fig. S2). We obtained 173 differentially expressed genes that could be used for prognosis by univariate Cox proportional hazards regression. Genes were subjected to the lymphocyte activation function, which was obtained from Gene Oncology (GO). This resulted in 12 lymphocyte activation-associated genes that were used for subsequent analysis. Finally, we obtained four genes (TCF3, IGFBP2, TYRO3 and NOD2) that were significantly related to overall survival using a stepwise regression algorithm. We used the lymphocyte activation-associated gene signature to calculate each patient’s risk score based the four genes’ expression levels weighted by regression coefficients in the univariate Cox proportional hazards regression analysis (Fig. S1). Positive coefficients of IGFBP2, TYRO3, and NOD2 (0.2992, 0.2421068 and 0.2155469) suggested that their expression was associated with poor prognosis in GBM patients, while the negative coefficient of TCF3 (−0.6246210) indicated it was associated with better survival.

In the training dataset, GBM patients were divided into high-risk (n = 262) and low-risk (n = 263) groups according to the median risk score (2.752) of the lymphocyte activation-associated gene signature. This score was used as the cutoff. Patients within high-risk group had significantly shorter overall survival time than those within low-risk group (median survival 15.3 months vs 12.6 months, HR = 1.55, 95% CI [1.28–1.87], log-rank test P < 0.001, Fig. 1A). The gene expression and survival time distributions suggested that patients in the high-risk group had higher IGFBP2 expression (Fig. 1B). The low-risk patients tended to have a longer survival time (Fig. 1C).

Prognostic value of the lymphocyte activation-associated gene signature in validation sets

Next, we validated the prognostic value of the lymphocyte activation-associated gene signature in another two external validation datasets, GSE13041 (n = 191) and GSE16011 (n = 150). Each patient’s risk score was calculated with the same formula as that used in the training set. The patients were divided into the high- and low-risk groups according to the median cutoff determined for each dataset. In GSE16011, survival time of patients in the high-risk group were significantly shorter than those in the low-risk group (median survival 7.1 months vs 10.3 months, HR = 1.57, 95% CI [1.13–2.17], log-rank test P = 0.069, Fig. 2A). IGFBP2 and TYRO3 had a relatively higher expression in the high-risk group (Fig. 2B). We found the low-risk patients presented a trend of longer survival time in GSE16011 (Fig. 2C). Moreover, the signature was found to predict the survival of patients with GBM in GSE13041 (median survival of high-risk patients 12.0 months vs median survival of low-risk patients 14.5 months, HR = 1.42, 95% CI [1.05–1.92], log-rank test P = 0.02, Figs. 3A–3C) with statistical significance.

The signature independently predicted overall survival of GBM patients

We used univariate and multivariable Cox regression analysis after adjusting for other clinicopathologic factors to assess whether the lymphocyte activation-associated gene signature was an independent prognostic biomarker in GBM. The covariables included age, sex, IDH status, subtypes, and the prognostic signature. We found that the signature (HR = 1.46, 95% CI [1.16–1.83], P = 0.001), age (HR = 1.03, 95% CI [1.02–1.04], P < 0.001), sex (HR = 1.30, 95% CI [1.04–1.63], P = 0.02) and proneural types (HR = 1.64, 95% CI [1.18–2.27], P = 0.003) independently predicted a worse OS for GBM patients in the training dataset (Table 3). In the validation dataset GSE16011, we also found that the lymphocyte activation-associated gene signature had independently prognostic value (HR = 1.36, 95% CI [1.35–1.98], P = 0.029, Table 4). This signature independently predicted the survival of GBM patients, with marginal significance (HR = 1.33, 95% CI [0.98–1.81], P = 0.06, Table S1). These findings revealed that the prognostic ability of the lymphocyte activation-associated gene signature was independent of the clinicopathological factors for OS in GBM.

The differences in lymphocyte activity between the high- and low-risk groups of patients

The GBM patients were divided into low- and high-risk groups based on their lymphocyte activation-associated gene signature. We then compared the lymphocyte activity in these two groups. Lymphocytes include natural killer cells, T cells and B cells. In TCGA, five types of lymphocytes showed significant differences in activity between high- and low-risk patients by Wilcoxon ranked sum test (P < 0.05). These lymphocytes were memory B cells, naive B cells, naive CD4 T cells, follicular helper T cells and regulatory T cells (Fig. 4A and Fig. S3). And the lymphocyte activity of high-risk patients was significantly lower than in low-risk patients. Similar phenomena were observed in the GSE13041 dataset (Fig. 4B). However, activated CD4 memory T cells, resting CD4 memory T cells, CD8 T cells and gamma delta T cells showed significant differences between the two groups. There were three types of cells that showed different activity between low- and high-risk patients in the GSE16011 dataset (P < 0.05, Fig. S4). These findings suggested that lymphocyte activity was significantly different between the high- and low-risk GBM patients.

Stratification analysis of the lymphocyte activation-associated gene signature

The IDH mutation is one of the most critical genomic alterations in GBM. The IDH mutation means a somatic mutation in IDH1 in our study. IDH-wildtype GBM patients had a shorter survival time than those with the IDH-mutation. The lymphocyte activation-associated gene signature significantly predicted the overall survival of IDH-wildtype GBM patients by log-rank test (P = 0.0043) but it did not predict the survival of IDH-mutation patients (Fig. 5A). The IDH-wildtype patients were divided into high-risk and low-risk groups using the same cutoff in the training dataset. The high-risk group of IDH-wildtype patients had significantly shorter OS than those in the low-risk group (median survival 12.9 months vs 14.9 months, HR = 1.39, 95% CI [1.11–1.74], P = 0.0043, Fig. 5A). The same phenomenon was also found in the GSE16011 dataset. The high-risk patients showed worse OS than the low-risk patients (median survival 7.1 months vs 10.7 months, HR = 1.95, 95% CI [1.25–3.02], P = 0.0025) among IDH-wildtype patients (Fig. 5B).

We next explored whether the lymphocyte activation-associated gene signature was effective for GBM patients within four transcriptome subtypes (proneural, neural, classical, and mesenchymal) established by Verhaak et al. using Kaplan–Meier survival analysis. The signature significantly predicted the overall survival of patients with the mesenchymal and proneural subtypes (log-rank test P value = 0.013 for mesenchymal subtypes and P value = 0.041 for proneural subtypes, Fig. 6). The patients in the high-risk groups had a shorter survival time than those in the low-risk groups (median survival 11.8 months vs 15.3 months for mesenchymal subtypes and 9 months vs 13.2 months for proneural subtypes, respectively). In addition, the signature did not have the prognostic ability for patients with neural and classical subtypes (Fig. S5).

We also explored whether radiotherapy was an effective treatment for GBM patients. We observed that high risk score significantly predicted a poor OS for GBM patients undergoing radiotherapy in the GSE16011 (median survival 8.76 months vs 14.64 months, HR = 1.86, 95% CI [1.27–2.72], P = 0.0011) and GSE13041 datasets (median survival 6.03 months vs 12.53 months, HR = 3.18, 95% CI [1.13–8.94], P = 0.022, Fig. 7).

Comparison of prognostic power between clinical factors and the signature

To evaluate the prognostic performance of the prognostic signature, we performed C-index calculation of clinical factors (age and sex) and the signature in three datasets (Harrell Jr et al, 1982; Pencina & D’Agostino, 2004). We constructed another prognostic model by combining our signature with clinical factors. A higher C-index indicated a superior prognostic value of signature (Wu, Yuan & Liu, 2019; Li, Cui & Li, 2017). There were significant differences found between age and the signature (P < 0.05) and between age and sex. The clinical factors for age in combination with the signature showed a higher C-index (0.646 for TCGA, 0.665 for GSE16011 and 0.600 for GSE13041) than age or histological grade alone (0.645 for TCGA, 0.659 for GSE16011 and 0.595 for GSE13041, Table S2), which were statistically significant (P < 0.05). These results suggested that this lymphocyte activation-associated gene signature could add complementary value to known clinical factors.