A novel immune-related gene prognostic signature combining immune cell infiltration and immune checkpoint for glioblastoma patients

Introduction

Cancer is a highly heterogeneous disease characterized by complex interactions between malignant cells and the tumor microenvironment (TME). The dynamic interaction of cancer cells with the microenvironment comprising cellular components (such as stromal cells and immune cells) and non-cellular components [extracellular matrix (ECM)] is essential to stimulate tumorigenesis, therapeutic response, as well as tumor progression and metastasis (1). Immune cells within the TME critically take part in the regulation of tumor growth and progression, therapeutic outcomes, and patient prognosis (2).

Gliomas represent a heterogeneous group of universally lethal brain tumors associated with dismal 5-year overall survival (OS) rates of less than 35% (3). It has long been suggested that the glial stem cell and/or progenitor cell population may be the origin of this disease and provide the underlying heterogeneity. On the basis of the histological phenotype, gliomas can be historically classified into astrocytomas, oligodendrogliomas, and ependymomas (4), and assigned grades I–IV by World Health Organization (WHO). Glioblastoma (GBM), a grade IV astrocytoma, is considered the most malignant intracranial tumor with a 5-year OS rate of less than 10% (5). GBM is notorious for rapid cell proliferation, angiogenesis, and therapy resistant glioblastoma stem-like cells (GSCs), leading to disease relapse and poor prognosis (6). The current standard therapy consists of surgery followed by radiotherapy and chemotherapy, but most cases eventually demonstrate chemoresistance to temozolomide (TMZ). Hence, identifying prognostic biomarkers associated with the biological heterogeneity of GBM is crucial for improving patient prognosis.

Immuno-oncology has been heralded as a significant breakthrough in cancer treatment. Emerging reports have revealed the existence of the lymphatic system within the central nervous system (CNS), pinpointing the vital role of immune regulation in treating brain tumors. Immunotherapy holds promise for a variety of primary cancers such as breast cancer and hepatocellular carcinoma (7). However, there is a persistent challenge for the application of immunotherapy to GBM due to the redundant mechanisms of tumor-mediated immunosuppression (8). Myeloid cells are significant components of the innate immune system within the GBM TME (9). While checkpoint blockade immunotherapy has been demonstrated to induce robust intratumoral T-cell infiltration in patients with brain metastases, these immunotherapeutic approaches fail to exert effective roles in GBM (10). Recent study has highlighted the development and validation of inflammatory response-related prognostic models and immune infiltration analyses in GBM, which underscores the importance of understanding immune responses in the TME and demonstrates the potential for immune-related biomarkers to improve prognosis and therapeutic strategies (11). Therefore, a comprehensive analysis of the relationship between immune-related genes and OS can offer new insights for the treatment and prognosis of GBM.

Cancer immunotherapy, as a novel therapeutic approach following surgery, chemotherapy, radiotherapy, and targeted therapy, has revolutionized the treatment paradigm for malignant tumors (12). The advancements in genome sequencing and bioinformatics have facilitated high-throughput analysis and interpretation of complex disease-related data, making them ideal for quantifying tumor-infiltrating immune cells in various cancers (13). Single-sample gene set enrichment analysis (ssGSEA) is an extension of gene set enrichment analysis (GSEA) and calculates an immune enrichment score for each sample-gene set pair. This method transforms the gene expression profile of a single sample into a gene-set enrichment profile, allowing researchers to characterize tumor-infiltrating immune cells within the TME using immunohistochemistry and flow cytometry (1).

In this study, we focused on constructing a prognostic signature using immune-related genes and assigned GBM patients into two clusters based on the results of ssGSEA. We identified 312 differentially expressed immune-related genes, of which 28 were correlated with the prognosis of GBM patients. Through least absolute shrinkage and selection operator (LASSO) regression analysis, 10 genes were screened out for model construction eventually. Correlation analysis between the signature and immune infiltrates, GBM-related genes, and tumor mutation burden (TMB) was performed. Finally, a nomogram was successfully established for prognosis prediction. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-562/rc).

Methods

Construction and validation of immune clustering

RNA sequencing data and relevant clinical information of GBM patients were obtained from The Cancer Genome Atlas (TCGA). By virtue of the gene set variation analysis (GSVA) package in R, ssGSEA was performed based on the expression levels of 29 immune cell subtypes. The GBM samples were grouped into the high and low immune infiltration clusters utilizing the “hclust” package in R. The hierarchical relationship and distribution of GBM samples were visualized using the “sparcl” and “Rtsne” packages. To validate the clustering results, we conducted differential analysis of immune infiltration between the two clusters. The “estimate” package in R was employed to evaluate the Immune Score, Stromal Score, ESTIMATE score, and Tumor Purity of GBM samples based on the immune infiltration scores of 29 immune cell subtypes. The “CIBERSORT” package was utilized to compare the expression levels of the human leukocyte antigen (HLA) family and 22 immune cell subtypes between the high and low immune infiltration clusters. Violin plots, heatmaps, and boxplots were created using the “reshape2”, “ggpubr”, and “pheatmap” packages to illustrate the findings. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

ssGSEA

ssGSEA between the two clusters was conducted using the GSEA software (v4.1.0). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) analysis were performed to compare the molecular mechanisms and cellular functions between the two clusters. The threshold for significance was set as false discovery rate (FDR) <0.001 for KEGG analysis and FDR <0.0001 for GO analysis. Bubble charts depicting the results were generated using the “reshape” and “ggplot2” packages.

Identification of differentially expressed immune-related genes between two clusters

Analysis of differentially expressed immune-related genes between the high and low immune clusters was carried out with a threshold of |log fold change (FC)| >0.585 and FDR <0.05. A list of immune-related genes was obtained from the Immunology Database and Analysis Portal (ImmPort, https://www.immport.org/), and Venn analysis was performed to illustrate the differentially expressed immune-related genes specific to the two clusters. The “ggplot2”, “venn”, and “pheatmap” packages in R were employed for visualization.

Correlation analysis between prognostic immune-related genes (PIGs) and transcription factors (TFs)

The gene expression files and clinical data of GSE83300 (platform GPL6480-9577) were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database (14). Univariate Cox regression analysis of the clinical data from TCGA was conducted to identify the PIGs correlated with the OS of GBM patients. A significance threshold of P<0.01 was applied. Co-expression analysis between PIGs and TFs downloaded from the Cistrome platform (http://www.cistrome.org/) was performed with thresholds of |cor| >0.4 and FDR <0.001. Sankey diagram was created with the “ggalluvial”, “ggplot2”, and “dplyr” packages in R. Additionally, a protein-protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, https://string-db.org/) online platform. A confidence score threshold of 0.9 (highest) was applied, and the disconnected nodes in the network were hidden.

Construction and validation of the immune-related prognostic signature

Using the identified PIGs, LASSO regression analysis was performed to mitigate overfitting during model construction. The prognostic model assigned a risk score to each patient based on the formula: Riskscore=∑i=1nCoePIGi×ExpPIGi, where Coe represents the coefficient of the PIG obtained from LASSO analysis and Exp represents the expression level of the PIG. The patients from TCGA and GSE83300 were divided into the high-risk and low-risk groups based on the median risk score, with TCGA patients serving as the training cohort and GSE83300 patients as the test cohort. Survival analysis with Kaplan-Meier (K-M) curves was performed in both cohorts, and time-dependent receiver operating characteristic (ROC) curves were drawn to validate the efficacy of the prognostic model. The R packages “survival”, “survminer”, “timeROC”, and “rms” were utilized for survival analysis. Furthermore, univariate and multivariate Cox regression analyses were conducted to determine whether the risk score was independent of other clinical features (such as age and gender) in predicting the prognosis for GBM patients. A significance level of P<0.05 was considered statistically significant.

Correlation analysis between genes for model construction and immune infiltrates

The correlation between PIGs in the model and immune cell infiltration was analyzed using the “CIBERSORT” algorithm, with a significance threshold of P<0.05. The R packages “tidyverse”, “ggplot2”, and “ggExtra” were used to visualize the results.

Correlation analysis between prognostic model and GBM-related genes

Differential expressed gene (DEG) analysis and correlation analysis between GBM-related genes and prognostic model were performed. Caveolae associated protein 1 (CAVIN1), isocitrate dehydrogenase (NADP (+)) 2 (IDH2), programmed cell death 1 (PDCD1), and cytotoxic T-lymphocyte associated protein 4 (CTLA4) were identified. The R packages “ggplot2”, “ggpubr” and “ggExtra” were used to plot the boxplots and correlation curves.

TMB analysis

TMB data of GBM patients were obtained from TCGA, and VarScan2 was used for data processing and further analysis. Differential analysis of TMB between the high-risk and low-risk groups was performed. Survival analysis was conducted to investigate the effect of TMB on OS of GBM patients. The patients were divided into the high-TMB and low-TMB groups according to the cutoff value at 1 mutation per megabase (mut/MB). The prognosis of patients in the two groups (high-TMB and low-TMB) and two risk groups (high-risk and low-risk) was analyzed, and K-M curves were applied to demonstrate the results. A value of P<0.05 was considered statistically significant. The R packages “survival” and “survminer” were involved in the procedure.

Construction and validation of nomogram

A practical nomogram was developed using the “rms” and “survival” packages in R, incorporating age, gender, and risk score obtained from the previous analysis as variables. The nomogram provides a graphical representation of the predictive model and allows for the estimation of individualized patient outcomes. Additionally, ROC curves were generated to assess the sensitivity and specificity of the nomogram, providing a measure of its performance in predicting outcomes.

Statistical analysis

The original data were summarized and processed using the Strawberry Perl v5.0.1. The Benjamini-Hochberg method was applied in differential analysis to control the FDR. Survival analysis for GBM patients was conducted using the K-M method. Statistical analyses were performed using the R software (version 4.1.1). A significance level of P<0.05 was indicative of a significant difference.

Results

Construction and validation of immune clustering

RNA sequencing data of 39,741 mRNAs from 169 GBM samples and corresponding clinical data of 599 GBM patients were obtained from TCGA database. According to the immune enrichment scores in ssGSEA, the GBM samples were divided into the high (Immune-H, n=99) and low (Immune-L, n=70) immune infiltration clusters. To verify the validity of the clustering, the Immune Score, Stromal Score, ESTIMATE Score, and Tumor Purity were calculated with the “estimate” algorithm. Compared with the Immune-H cluster, the Immune-L cluster showed lower Immune Score, Stromal Score, and ESTIMATE Score, but higher Tumor Purity (Figure 1A,1B). Furthermore, “CIBERSORT” algorithm revealed that the Immune-H cluster had significantly higher expression levels in all 24 subtypes of the HLA family than the Immune-L cluster (Figure 1C). Higher proportions of T cells gamma delta and monocytes, and lower proportions of natural killer (NK) cells resting and macrophages M0 were found in the Immune-H cluster compared with the Immune-L cluster (Figure 1D).

Figure 1 Construction and validation of immune clustering. (A,B) Heatmap and violin plots showing the Immune Score, Stroma Score, ESTIMATE Score, and Tumor Purity of two clusters according to the “estimate” algorithm. (C) Boxplots displaying the differentially expressed genes of the HLA family between two clusters. (D) Boxplots showing the proportions of 22 immune cell subtypes of two clusters according to the “CIBERSORT” algorithm. *, P<0.05; **, P<0.01; ***, P<0.001. MHC, major histocompatibility complex; pDCs, plasmacytoid dendritic cells; IFN, interferon; aDCs, activated dendritic cells; APC, antigen-presenting cells; TIL, tumor infiltrating lymphocytes; Treg, regulatory T cells; CCR, C-C chemokine receptor; HLA, human leukocyte antigen; iDCs, immature dendritic cells; DCs, dendritic cells; Tfh, T follicular helper cells; NK, natural killer; TME, tumor microenvironment; Immunity-L, low immune infiltration; Immunity-H, high immune infiltration.

ssGSEA

ssGSEA was performed within the high and low immune infiltration clusters. GO analysis revealed that the differences between the two clusters were mainly reflected in positive regulation of cell killing, negative regulation of leukocyte-mediated immunity, macrophage migration, etc. (Figure 2A). KEGG analysis demonstrated that primary immunodeficiency, intestinal immune network for immunoglobulin A (IgA) production, etc. were involved in distinguishing the two clusters (Figure 2B).

Figure 2 Bubble charts showing the results of Gene Set Enrichment Analysis. (A) Top 16 results of GO analysis. (B) Top 11 results of KEGG analysis. FDR, false discovery rate; NES, normalized enrichment score; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of differentially expressed immune-related genes between two clusters

With a threshold of |logFC| >0.585 and FDR <0.05, a total of 1,855 DEGs between the two clusters were obtained, including 859 genes with lower expression and 996 genes with higher expression in the Immune-H cluster (Figure 3A). However, no significant difference was found in the other 37,885 genes. A list of 1,793 immune-related genes was obtained from ImmPort. Venn analysis of the two sets of genes identified 312 immune-related DEGs (Figure 3B), including 31 DEGs with higher expression in the Immune-L cluster and 281 DEGs with higher expression in the Immune-H cluster (Figure 3C).

Figure 3 Identification of differentially expressed immune-related genes. (A) Volcano plot showing 1,855 DEGs, including 859 genes with lower expression (green) and 996 genes with higher expression (red) in the Immune-H cluster. (B) Venn diagram showing the intersection of DEGs and immune-related genes (n=312). (C) Heatmap showing 31 DEGs with higher expression in the Immune-L cluster and 281 DEGs with lower expression in the Immune-H cluster. FDR, false discovery rate; FC, fold change; DEGs, differential expressed genes; Immunity-L, low immune infiltration; Immunity-H, high immune infiltration.

Establishment of functional interaction networks of PIGs and TFs

The gene expression profiles of 50 GBM patients were obtained from GSE83300, covering 19,593 genes. Through co-expression analysis, all 312 differentially expressed immune-related genes were found to be shared by both TCGA and GSE83300 datasets. Univariate Cox regression analysis of the clinical data from both datasets showed that a total of 28 genes were correlated with the OS of patients (Figure 4A). From the Cistrome platform, a total of 317 TFs were obtained, among which 19 TFs, including BCL3, BATF, CEBPB, etc., were found to be related to the 28 PIGs (Figure 4B). Furthermore, the PPI network between the 19 TFs and 28 PIGs was plotted using STRING, with FOS, JUN, STAT1, IRF1, etc. being central to the network (Figure 4C).

Figure 4 Correlation between PIGs and TFs. (A) Forest map demonstrating 28 genes related to prognosis. (B) Sankey diagram showing the correlation between 28 PIGs and 19 TFs. (C) A protein-protein interaction network consisting of PIGs and TFs. PIGs, prognostic immune-related genes; TFs, transcription factors.

Construction and validation of the immune-related prognostic signature

LASSO regression analysis identified 10 out of 28 PIGs for model construction, including CLCF1, PTX3, TNFRSF14, SDC2, VGF, AREG, PLAUR, GRN, AQP9, and IGLV6-57 (Figure 5A,5B). The coefficients of the 10 PIGs are shown in Table 1. A total of 594 patients from the TCGA dataset were assigned as the training cohort, while 50 patients from the GSE83300 dataset were used as the test cohort. Based on the expression levels and coefficients of the 10 genes, each patient was assigned a risk score and then allocated to the high-risk group or the low-risk group. Survival analysis revealed that patients in the low-risk group had a better prognosis compared to those in the high-risk group (P<0.001 in the training cohort, P=0.007 in the test cohort; Figure 5C,5D). To evaluate the sensitivity and specificity of the prognostic model, ROC analysis was conducted. The area under the curve (AUC) values at 1, 3, and 5 years were 0.750, 0.770, and 0.819 in the training cohort, and 0.637, 0.721, and 0.659 in the test cohort, respectively (Figure 5E,5F). Additionally, Figure 5G,5H demonstrate a negative correlation between OS and risk score in both the training and test cohorts (R=−0.38 in the training group, R=−0.24 in the test group). The risk score distribution is shown in Figure 5I,5J, while Figure 5K,5L display the survival status of each GBM patient in the training and test cohorts. Univariate and multivariate Cox regression analyses indicated that the risk score was a prognostic factor independent of other clinical characteristics (Figure 5M,5N). Although the P value for age was less than 0.001, indicating statistical significance, the hazard ratio (HR) was 1.026. This suggests that while age has a statistically significant association with the outcome, its effect size is minimal, and thus, it does not have a substantial impact on the overall prognosis.

Figure 5 Construction and validation of an immune-related prognostic model. (A,B) Plots showing the LASSO coefficient profiles and partial likelihood deviance for 28 PIGs, and 10 genes were selected for model construction. (C,D) Patients in the low-risk group showed better OS than those in the high-risk group in both the training cohort (P<0.001) and the test cohort (P=0.007). (E,F) ROC curves of the prognostic model at 1, 3, and 5 years in the training cohort (AUC =0.750, 0.770, 0.819) and the test cohort (AUC =0.637, 0.721, 0.659). (G,H) The OS of GBM patients was negatively correlated with the risk score in the training cohort (R=−0.38) and the test cohort (R=−0.24). (I-L) Distribution of the risk score and survival status of patients in the training cohort (I,J) and the test cohort (K,L). (M,N) Forest plots of univariate and multivariate Cox analyses showing that the risk score was an independent predictor of OS for GBM patients, independent of age and gender (P<0.001). AUC, area under the curve; OS, overall survival; LASSO, least absolute shrinkage and selection operator; PIGs, prognostic immune-related genes; ROC, receiver operating characteristic; GBM, glioblastoma.

Correlation between PIGs in the model and immune infiltrates

To examine the correlation between the 10 genes in the prognostic model and immune function, we conducted correlation analysis using the “CIBERSORT” algorithm. The results revealed that IGLV6-57 had a strong correlation with immune cells. Particularly, CD4 memory activated T cells showed a positive correlation with most of the genes, and neutrophils exhibited a particularly strong positive correlation with AQP9 (Figure 6).

Figure 6 Correlation between the genes in the prognostic model and immune infiltrates. NK, natural killer.

Correlation between GBM-related genes and prognostic signature

Through differential expression analysis and correlation analysis, we observed that the high-risk group presented higher expression levels of CAVIN1 (Figure 7A) and lower expression levels of IDH2 (Figure 7B). CAVIN1 was positively correlated with the risk score (Figure 7C), while IDH2 showed a negative correlation with the risk score (Figure 7D). Similarly, we found higher expression levels of PDCD1 (Figure 7E) and CTLA4 (Figure 7F), two immune checkpoint-related genes, in the high-risk group. PDCD1 and CTLA4 were positively correlated with the risk score (Figure 7G,7H).

Figure 7 Correlation between the prognostic model and GBM-related immune genes. (A-D) CAVIN1 was highly expressed in the high-risk group and positively correlated with the risk score. IDH2 was poorly expressed in the high-risk group and negatively correlated with the risk score. (E-H) PDCD1 was highly expressed in the high-risk group and positively correlated with the risk score. CTLA4 was highly expressed in the high-risk group and positively correlated with the risk score. CAVIN1, caveolae associated protein 1; IDH2, isocitrate dehydrogenase [NADP (+)] 2; PDCD1, programmed cell death 1; CTLA4, cytotoxic T-lymphocyte associated protein 4; GBM, glioblastoma.

TMB analysis

TMB data of GBM patients were processed with VarScan2 and downloaded from TCGA. Differential analysis revealed that the high-risk group had lower TMB compared to the low-risk group (P=0.01; Figure 8A). With a cutoff value of 1 mut/Mb, patients were divided into the high-TMB and low-TMB groups. Survival analysis showed that patients in the high-TMB group tended to have better OS (P=0.06, Figure 8B) compared with those in the low-TMB group. Furthermore, patients in the high-TMB and low-risk groups exhibited significantly better prognoses compared to those in the low-TMB and high-risk groups (P<0.001, Figure 8C).

Figure 8 TMB analysis. (A) Boxplots showing higher TMB levels in the low-risk group (P=0.01). (B) Kaplan-Meier curves indicating that patients in the high-TMB group had better OS compared to those in the low-TMB group (P=0.06). (C) Patients in the high-TMB and low-risk groups showed better prognosis than those in the low-TMB and high-risk groups (P<0.001). H, high; L, low; TMB, tumor mutation burden; OS, overall survival.

Construction and validation of nomogram

To provide a prediction model with enhanced practicality and validity, a nomogram incorporating age, gender, and risk score was constructed for clinical application (Figure 9A). The AUC values at 1, 3, and 5 years for the nomogram were 0.767, 0.870, and 0.841, respectively, indicating better validity compared to the prognostic signature alone (Figure 9B).

Figure 9 Construction and validation of nomogram. (A) A nomogram incorporating age, gender, and risk score was established. (B) ROC curves of the nomogram for 1-, 3-, and 5-year OS prediction (AUC =0.767, 0.870, and 0.841). ***, P<0.001. AUC, area under the curve; OS, overall survival; ROC, receiver operating characteristic.

Discussion

GBM exhibit extensive inter- and intra-tumoral heterogeneity, complicating the disease and therapeutic strategies (15). The immune infiltration characteristics in the TME have been identified as potential sources of prognostic biomarkers. However, inconsistent results and limited translational impact are observed due to the heterogeneity in patient populations, treatment regimens, and technology platforms (2). In recent years, next-generation sequencing data available through public databases such as TCGA and GEO have enabled risk stratification based on the highly heterogeneous nature of the tumor. In this study, we constructed a risk model consisting of 10 immune-related genes using bioinformatics methods. These immune-related genes were found to be strongly correlated with the prognosis and immune infiltration of GBM patients (Figure 6). GBM samples from the low-risk group exhibited higher expression levels of genes associated with adaptive immune response, yielding better prognoses. The immune-related model not only improved the prognostic predictive ability but also provided insights into the immune mechanisms underlying tumorigenesis. Additionally, we observed higher expression levels of CAVIN1 and PDCD1, but lower levels of IDH2 and CTLA4 in tumor samples from the high-risk group (Figure 7), suggesting that these genes might have potential correlations with the risk score. By incorporating clinical parameters including age, gender, and risk score, we developed a nomogram for predicting OS of GBM patients (Figure 9). Both the risk prediction model and nomogram showed good accuracy, reliability, and sensitivity in view of the ROC curves.

To explore the differences between the high and low immune infiltration clusters as well as the potential mechanisms, we used ssGSEA to perform GO and KEGG analyses on the genes in the target database. The existing studies have revealed the vital roles of macrophage migration (16), negative regulation of immune system process (17), and T cell proliferation (18) in the malignant process of GBM, which is confirmed by our results. We observed that 7 of the 10 immune-related genes (AQP9, AREG, CLCF1, GRN, PLAUR, PTX3, and SDC2) were highly expressed in the GBM group, while two of these immune-related genes (TNFRSF14 and VGF) were poorly expressed. However, IGLV6-57 was not clearly expressed in GBM. Most of the 10 immune-related genes have been associated with the prognosis of different types of cancers. AQP9 is not only a water channel for transporting glycerol, mannitol, and urea, but also a peroxiporin that efficiently transports H₂O₂ (19). Under normal conditions, AQP9 is weakly expressed in the brain, while widespread AQP9 expression has been reported in human GBM (20). In GBM, most cells throughout the tumor show strong anti-AQP9 immunoreactivity on the whole cell surface.

GBM cells often undergo metabolic reprogramming to support rapid growth and survival, including aerobic glycolysis (a Warburg-like effect) and other metabolic dysfunctions such as amino acid metabolism and lipid metabolism (21). Immune-related gene-classified GBM may show distinct metabolic profiles, particularly in the context of immune cell infiltration and activity. For example, higher glycolytic activity in tumor cells can lead to a hypoxic environment, influencing immune cell function and recruitment (22). The hypoxic environment in GBM can lead to the accumulation of metabolites like lactate, which can suppress immune responses by affecting the metabolism and function of immune cells (23). Certain metabolic pathways, such as those involving tryptophan and arginine, can be modulated to affect immune cell activity within the TME (24). As a known gene involved in metabolic pathways, AQP9 can facilitate glycerol and lactic acid clearance, which is critical for maintaining the metabolic balance in GBM cells. The enhanced immunoreactivity against AQP9 in GBM cells is thought to counteract GBM-associated lactic acidosis by clearing glycerol and lactic acid from the extracellular space. In addition, the augmented AQP9 immunoreactivity may be associated with the energy metabolism of glioma and/or surrounding neuronal cells (25). The ability of tumor stem cells to migrate around normal tissues can be partially attributed to the high expression of AQP9 aquaporin, as some cells acquire a more “differentiated” phenotype (26). Our results unveiled a highly positive correlation between AQP9 and activated memory CD4 T cells and neutrophils. It is suggested that AQP9 may also be a candidate for targeted immunotherapy, since at least its cell surface should be accessible to antibodies. In addition to AQP9, AREG, IGLV6-57, and PLAUR were also highly positively correlated with activated memory CD4 T cells (Figure 6). AREG, a member of the epidermal growth factor receptor (EGFR) family and a ligand of EGFR (27), is known as a bi-functional growth-regulatory glycoprotein that stimulates several human fibroblast cell lines while inhibiting some neuroblastoma cell lines (28). Alterations of AREG receptor in malignant astrocytomas are responsible for the activation of a complex network of pathways including Ras/MAPK, PI3K/AKT, PLCγ, and STAT (29). AREG has been demonstrated as a gene associated with glioma cell migration and invasion in several genome-wide studies. AREG gene promoter methylation is significantly correlated with poor prognosis of GBM. Moreover, AREG methylation is identified as an independent prognostic predictor in Cox multivariate analysis (30). The function of IGLV6-57 in GBM has rarely been reported, and we are the first to demonstrate the immune role of IGLV6-57 in GBM. IGLV6-57 was highly positively correlated with naive CD8 T cells and B cells, while negatively correlated with regulatory T cells (Tregs) and M0 macrophages (Figure 6). PLAUR encodes the urokinase receptor [urokinase-type plasminogen activator receptor (uPAR)], which promotes cell survival, migration, and resistance to targeted therapies in GBM cells. PLAUR gene expression adversely affects GBM patient survival by promoting mesenchymal gene expression profiles, allowing cell survival, and inducing stem-cell-like properties in a small population of GBM cells. A recent study has identified PLAUR as a valuable prognostic biomarker in glioma (31). CLCF1, a neurotrophic and B cell stimulator belonging to the interleukin-6 (IL-6) family (32), is found to induce the phosphorylation of signal transducers and transcriptional activators 3 (STAT3). Continuous activation of STAT3 in phosphatase and tensin homolog (PTEN)-deficient GBM induces cell proliferation, anti-apoptosis, maintenance of tumor stem cells, tumor invasion, angiogenesis, and immune evasion. Inhibition of CLCF1 to reduce STAT3 phosphorylation may be an effective therapeutic strategy for PTEN-mutated gliomas (33). As a protein coding gene, GRN plays an important role in tumorigenesis and functions as a prognostic marker of GBM. Researchers have demonstrated the association between increased GRN expression and poor prognosis in GBM (34). Our results also showed that GRN was positively correlated with activated memory CD4 T cells and activated dendritic cells, indicating the potential of GRN as a target of immunotherapy (Figure 6). PTX3 actively mediates the infiltration, migration, and anti-inflammatory polarization of macrophages, and also predicts immunotherapy responses in GBM (35). PTX3 expression positively correlates with glioma grade (36). The high expression of PTX3 in GBM, either caused by the tumor itself or by the neighboring cells, may contribute to the activation of the inflammatory program in the TME. In GBM, PTX3 immunoreactivity is stronger in cells closest to the vascular network or in palisade cells surrounding the necrotic area (37), reflecting tumor inflammation and malignancy. Our results showed that SDC2 was related to many immune factors (Figure 6). As a member of the syndecan family, SDC2 is constitutively expressed in normal brain tissues and GBM (38), while its expression in GBM specimens is extremely higher than that in normal brain specimens. However, there is a lack of research on the role of SDC2 in glioma immunity, and our model can arouse attention on this target. Tumour necrosis factor receptor superfamily 14 (TNFRSF14) is another essential immune checkpoint molecule critical in immune surveillance for cancer, and has been suggested as a tumor suppressor gene in GBM (39). TNFRSF14 is weakly expressed in GBM and contributes to T-cell immunity suppression. Our results confirmed the association between TNFRSF14 and dendritic cells resting, which may be a potential tumor-promoting mechanism (40). Intriguingly, higher expression of TNFRSF14 is also found to be associated with poorer OS and disease-free survival and lower Th1 response (41). All these findings support TNFRSF14 as a potential prognostic marker or molecular target in the clinical treatment of GBM. VGF is a neural growth factor-induced gene that participates in the regulation of energy homeostasis, neurogenesis, synapse formation, and psychiatric disorders (42). Differentiated tumor cells synergistically interact with stem-like tumor cells to promote GBM hierarchy and tumor growth through a paracrine feedback loop of VGF-mediated neurotrophic signaling (43). Our results suggested the role of VGF in GBM immunotherapy, as evidenced by its regulation on monocytes and M0 macrophages (Figure 6), which has not been reported yet.

High TMB can lead to immune checkpoint blockade (ICB) responses in several cancers, with glioma being the notable exception (44). Common gene mutations in GBM include TP53, PTEN, EGFR, and IDH1/2. These mutations can influence the immune landscape by affecting the expression of immune checkpoints and the presentation of neoantigens (45,46). Immune-related gene-classified GBM may show specific mutation signatures that correlate with immune evasion mechanisms, such as upregulation of programmed cell death 1 ligand 1 (PD-L1) or alterations in antigen presentation machinery (47,48). High mutational burden can lead to the generation of neoantigens, which can make tumors more immunogenic and potentially responsive to immunotherapy (49). Several clinical immunotherapy trials have reported long-term survival benefits in 10–20% of patients with recurrent GBM. Specifically, a lower TMB is associated with a longer survival of patients with recurrent GBM after recombinant poliovirus therapy or ICB. However, no association between TMB and survival is observed in the cohort of patients with newly diagnosed or recurrent GBM patients who do not receive immunotherapy. In the cohort of recurrent but not newly diagnosed GBM, there is an inverse relationship between TMB and enrichment of inflammatory gene signatures, suggesting that the relationship between TMB and inflammation inherent in the tumor evolves at the time of recurrence (50). In our analysis, we found that patients with high-TMB survived longer than those with low-TMB. Moreover, the TMB of patients in the low-risk group was higher than that in the high-risk group, which also resulted in the longer survival time of patients with high-TMB and low-risk. Future evaluation of the correlation between TMB and prognosis in different GBM subgroups is expected to validate our findings and guide the application of immunotherapy.

The immune profile of GBM includes the presence of tumor-infiltrating lymphocytes (TILs), macrophages, microglia, and other immune cells. The composition and functional state of these cells can significantly impact tumor cell behaviors and therapeutic responses (51). Immune-related gene classification can reveal differences in the abundance and activation state of immune cells. For instance, a high presence of Tregs or myeloid-derived suppressor cells (MDSCs) can indicate an immunosuppressive environment (52,53). The expression of cytokines and chemokines in GBM can modulate immune cell recruitment and function. For example, high levels of IL-10 and TGF-β are indicative of immunosuppression (54). Immune-related gene expression profiles contribute to identifying cytokine signatures that can predict the immune microenvironment and therapeutic response.

The initiation of T-cell antigen receptor (TCR) signaling is a key step that leads to the activation of T cells and the choreography of adaptive immune responses (55). The TME can utilize tyrosine kinase signaling as a molecular immune-rheostat to control immune responses by simultaneously altering T-cell activation and T-cell recruiting gene expression, especially the ability of the immune system to recognize and eradicate tumor cells (56). TNFRSF14 is upregulated following T cell priming and expressed on activated T cells, functioning as a co-regulatory molecule to mediate TCR signaling activation (57). SDC2, as a new regulator of v-domain Ig suppressor of T cell activation (VISTA) and monocyte binding, also participates in the TCR signaling process (58). GRN promotes immune evasion in pancreatic ductal carcinoma and serves as an efficient target to enhance cancer-specific antigen-mediated cytotoxicity for TCR (59). Accordingly, we speculate that GRN may also exert the same function in GBM, which needs to be proven by subsequent studies.

Immune-related genes can serve as biomarkers for predicting response to immunotherapies, such as checkpoint inhibitors [e.g., anti-programmed cell death 1 (PD-1)/PD-L1, anti-CTLA4] (60). For instance, patients with higher PD-L1 expression or specific immune gene signatures may benefit more from ICB therapy therapies (61). Exploring the expression profiles of immune-related genes can help in designing combination therapies that target both the tumor and its immune environment. For example, combining immune checkpoint inhibitors with therapies that modulate the TME (e.g., anti-angiogenic agents) can enhance therapeutic efficacy (62). Immune-related genes can also guide the use of therapies that modulate the immune system, such as vaccines, oncolytic viruses, or adoptive cell therapies [e.g., chimeric antigen receptor-T (CAR-T) cells] (63). Integrating immune-related gene signatures with other molecular and clinical data can promote the development of personalized treatment strategies, which enable patients to receive the most effective treatment based on their unique tumor and immune characteristics.

Finally, we designed a nomogram to predict the prognosis of GBM patients based on the risk score of 10 immune-related gene signature, gender, and age. From this nomogram we were able to predict the 1-, 3-, and 5-year survival of GBM. Consistently, the plots of calibration validated that the signature could accurately evaluate the survival of GBM patients. Despite the favorable predictive capability, the prognostic model can be further optimized by consolidating with other independent datasets or even improved by optimizing the LASSO results in the future.

There are some limitations in this study. The amount of data we downloaded from the public databases was relatively small, lacking sufficient sample size and other clinicopathological features, and moreover, publication bias and bulk effects could not be fully measured. Although our prognostic predictive features show favorable predictive power, objectively, the prediction model can be optimized by integrating with other independent datasets, and our prognostic features should also be validated in more GBM datasets. Due to limited conditions, the reliability of conclusions drawn from the model failed to be validated through our own clinical samples. Indeed, practical experiments are desirable to further verify the results based on public datasets. Last but not the least, the five parameters for constructing the nomogram may not be optimal due to the lack of other online parameters. Given the complex clinical characteristics of GBM, a better model and nomogram shall be built.

Conclusions

In conclusion, we developed a robust prognostic model for GBM in which the 10 immune-related genes were identified as independent prognostic predictors. Furthermore, a comprehensive nomogram was developed to reliably predict the OS of GBM patients. A better understanding of the relationship between immune-related genes, TMB, and immune infiltrates may provide insights into the underlying mechanisms of GBM tumorigenesis and guide individualized treatment. This study paves a novel way to predict the prognosis and survival of GBM and provides some promising targets for immunotherapy.

Acknowledgments

Funding: None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-562/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-562/prf

Conflicts of Interest: Both authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-562/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Zuo S, Wei M, Wang S, et al. Pan-Cancer Analysis of Immune Cell Infiltration Identifies a Prognostic Immune-Cell Characteristic Score (ICCS) in Lung Adenocarcinoma. Front Immunol 2020;11:1218. [Crossref] [PubMed]
2. Gentles AJ, Newman AM, Liu CL, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med 2015;21:938-45. [Crossref] [PubMed]
3. Lapointe S, Perry A, Butowski NA. Primary brain tumours in adults. Lancet 2018;392:432-46. [Crossref] [PubMed]
4. Weller M, Wick W, Aldape K, et al. Glioma. Nat Rev Dis Primers 2015;1:15017. [Crossref] [PubMed]
5. Campian JL, Ghosh S, Kapoor V, et al. Long-Acting Recombinant Human Interleukin-7, NT-I7, Increases Cytotoxic CD8 T Cells and Enhances Survival in Mouse Glioma Models. Clin Cancer Res 2022;28:1229-39. [Crossref] [PubMed]
6. Peng Z, Liu C, Wu M. New insights into long noncoding RNAs and their roles in glioma. Mol Cancer 2018;17:61. [Crossref] [PubMed]
7. Xiao B, Liu L, Li A, et al. Identification and Verification of Immune-Related Gene Prognostic Signature Based on ssGSEA for Osteosarcoma. Front Oncol 2020;10:607622. [Crossref] [PubMed]
8. Chen Q, Han B, Meng X, et al. Immunogenomic analysis reveals LGALS1 contributes to the immune heterogeneity and immunosuppression in glioma. Int J Cancer 2019;145:517-30. [Crossref] [PubMed]
9. Lee AH, Sun L, Mochizuki AY, et al. Neoadjuvant PD-1 blockade induces T cell and cDC1 activation but fails to overcome the immunosuppressive tumor associated macrophages in recurrent glioblastoma. Nat Commun 2021;12:6938. [Crossref] [PubMed]
10. Schaettler MO, Richters MM, Wang AZ, et al. Characterization of the Genomic and Immunologic Diversity of Malignant Brain Tumors through Multisector Analysis. Cancer Discov 2022;12:154-71. [Crossref] [PubMed]
11. Zhu W, Luo N, Li Q, et al. Development and validation of an inflammatory response-related prognostic model and immune infiltration analysis in glioblastoma. Ann Transl Med 2023;11:69. [Crossref] [PubMed]
12. Yuan F, Cai X, Zhu J, et al. A Novel Immune Classification for Predicting Immunotherapy Responsiveness in Patients With Adamantinomatous Craniopharyngioma. Front Neurol 2021;12:704130. [Crossref] [PubMed]
13. Newman AM, Alizadeh AA. High-throughput genomic profiling of tumor-infiltrating leukocytes. Curr Opin Immunol 2016;41:77-84. [Crossref] [PubMed]
14. Feng L, Qian H, Yu X, et al. Heterogeneity of tumor-infiltrating lymphocytes ascribed to local immune status rather than neoantigens by multi-omics analysis of glioblastoma multiforme. Sci Rep 2017;7:6968. [Crossref] [PubMed]
15. Stumpo V, Guida L, Bellomo J, et al. Hemodynamic Imaging in Cerebral Diffuse Glioma-Part B: Molecular Correlates, Treatment Effect Monitoring, Prognosis, and Future Directions. Cancers (Basel) 2022;14:1342. [Crossref] [PubMed]
16. Zhou Y, Meng X, He W, et al. USF1/CD90 signaling in maintaining glioblastoma stem cells and tumor-associated macrophages adhesion. Neuro Oncol 2022;24:1482-93. [Crossref] [PubMed]
17. Di Tacchio M, Macas J, Weissenberger J, et al. Tumor Vessel Normalization, Immunostimulatory Reprogramming, and Improved Survival in Glioblastoma with Combined Inhibition of PD-1, Angiopoietin-2, and VEGF. Cancer Immunol Res 2019;7:1910-27. [Crossref] [PubMed]
18. Diamant G, Simchony Goldman H, Gasri Plotnitsky L, et al. T Cells Retain Pivotal Antitumoral Functions under Tumor-Treating Electric Fields. J Immunol 2021;207:709-19. [Crossref] [PubMed]
19. Watanabe S, Moniaga CS, Nielsen S, et al. Aquaporin-9 facilitates membrane transport of hydrogen peroxide in mammalian cells. Biochem Biophys Res Commun 2016;471:191-7. [Crossref] [PubMed]
20. Jelen S, Parm Ulhøi B, Larsen A, et al. AQP9 expression in glioblastoma multiforme tumors is limited to a small population of astrocytic cells and CD15(+)/CalB(+) leukocytes. PLoS One 2013;8:e75764. [Crossref] [PubMed]
21. Garcia JH, Jain S, Aghi MK. Metabolic Drivers of Invasion in Glioblastoma. Front Cell Dev Biol 2021;9:683276. [Crossref] [PubMed]
22. Chen Z, Han F, Du Y, et al. Hypoxic microenvironment in cancer: molecular mechanisms and therapeutic interventions. Signal Transduct Target Ther 2023;8:70. [Crossref] [PubMed]
23. Park JH, Lee HK. The Role of Hypoxia in Brain Tumor Immune Responses. Brain Tumor Res Treat 2023;11:39-46. [Crossref] [PubMed]
24. Bader JE, Voss K, Rathmell JC. Targeting Metabolism to Improve the Tumor Microenvironment for Cancer Immunotherapy. Mol Cell 2020;78:1019-33. [Crossref] [PubMed]
25. Warth A, Mittelbronn M, Hülper P, et al. Expression of the water channel protein aquaporin-9 in malignant brain tumors. Appl Immunohistochem Mol Morphol 2007;15:193-8. [Crossref] [PubMed]
26. Fossdal G, Vik-Mo EO, Sandberg C, et al. Aqp 9 and brain tumour stem cells. ScientificWorldJournal 2012;2012:915176. [Crossref] [PubMed]
27. Yu Z, Yu Q, Xu H, et al. IL-17A Promotes Psoriasis-Associated Keratinocyte Proliferation through ACT1-Dependent Activation of YAP-AREG Axis. J Invest Dermatol 2022;142:2343-52. [Crossref] [PubMed]
28. Urbanavičiūtė R, Skauminas K, Skiriutė D. The Evaluation of AREG, MMP-2, CHI3L1, GFAP, and OPN Serum Combined Value in Astrocytic Glioma Patients' Diagnosis and Prognosis. Brain Sci 2020;10:872. [Crossref] [PubMed]
29. Steponaitis G, Kazlauskas A, Skiriute D, et al. Significance of Amphiregulin (AREG) for the Outcome of Low and High Grade Astrocytoma Patients. J Cancer 2019;10:1479-88. [Crossref] [PubMed]
30. Skiriutė D, Vaitkienė P, Ašmonienė V, et al. Promoter methylation of AREG, HOXA11, hMLH1, NDRG2, NPTX2 and Tes genes in glioblastoma. J Neurooncol 2013;113:441-9. [Crossref] [PubMed]
31. Wang LJ, Lv P, Lou Y. A Novel TAF-Related Signature Based on ECM Remodeling Genes Predicts Glioma Prognosis. Front Oncol 2022;12:862723. [Crossref] [PubMed]
32. Shen SH, Guo JF, Huang J, et al. Bromodomain-containing protein 4 activates cardiotrophin-like cytokine factor 1, an unfavorable prognostic biomarker, and promotes glioblastoma in vitro. Ann Transl Med 2022;10:475. [Crossref] [PubMed]
33. Zhang P, Meng X, Liu L, et al. Identification of the Prognostic Signatures of Glioma With Different PTEN Status. Front Oncol 2021;11:633357. [Crossref] [PubMed]
34. Vachher M, Arora K, Burman A, et al. NAMPT, GRN, and SERPINE1 signature as predictor of disease progression and survival in gliomas. J Cell Biochem 2020;121:3010-23. [Crossref] [PubMed]
35. Zhang H, Wang Y, Zhao Y, et al. PTX3 mediates the infiltration, migration, and inflammation-resolving-polarization of macrophages in glioblastoma. CNS Neurosci Ther 2022;28:1748-66. [Crossref] [PubMed]
36. Wesley UV, Sutton I, Clark PA, et al. Enhanced expression of pentraxin-3 in glioblastoma cells correlates with increased invasion and IL8-VEGF signaling axis. Brain Res 2022;1776:147752. [Crossref] [PubMed]
37. Locatelli M, Ferrero S, Martinelli Boneschi F, et al. The long pentraxin PTX3 as a correlate of cancer-related inflammation and prognosis of malignancy in gliomas. J Neuroimmunol 2013;260:99-106. [Crossref] [PubMed]
38. Watanabe A, Mabuchi T, Satoh E, et al. Expression of syndecans, a heparan sulfate proteoglycan, in malignant gliomas: participation of nuclear factor-kappaB in upregulation of syndecan-1 expression. J Neurooncol 2006;77:25-32. [Crossref] [PubMed]
39. Zhang Y, Cruickshanks N, Yuan F, et al. Targetable T-type Calcium Channels Drive Glioblastoma. Cancer Res 2017;77:3479-90. [Crossref] [PubMed]
40. Han MZ, Wang S, Zhao WB, et al. Immune checkpoint molecule herpes virus entry mediator is overexpressed and associated with poor prognosis in human glioblastoma. EBioMedicine 2019;43:159-70. [Crossref] [PubMed]
41. Lombardo SD, Bramanti A, Ciurleo R, et al. Profiling of inhibitory immune checkpoints in glioblastoma: Potential pathogenetic players. Oncol Lett 2020;20:332. [Crossref] [PubMed]
42. Muthukrishnan SD, Alvarado AG, Kornblum HI. Building Bonds: Cancer Stem Cells Depend on Their Progeny to Drive Tumor Progression. Cell Stem Cell 2018;22:473-4. [Crossref] [PubMed]
43. Wang X, Prager BC, Wu Q, et al. Reciprocal Signaling between Glioblastoma Stem Cells and Differentiated Tumor Cells Promotes Malignant Progression. Cell Stem Cell 2018;22:514-528.e5. [Crossref] [PubMed]
44. Samstein RM, Lee CH, Shoushtari AN, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet 2019;51:202-6. [Crossref] [PubMed]
45. Verdugo E, Puerto I, Medina MÁ. An update on the molecular biology of glioblastoma, with clinical implications and progress in its treatment. Cancer Commun (Lond) 2022;42:1083-111. [Crossref] [PubMed]
46. Decraene B, Vanmechelen M, Clement P, et al. Cellular and molecular features related to exceptional therapy response and extreme long-term survival in glioblastoma. Cancer Med 2023;12:11107-26. [Crossref] [PubMed]
47. Segura-Collar B, Hiller-Vallina S, de Dios O, et al. Advanced immunotherapies for glioblastoma: tumor neoantigen vaccines in combination with immunomodulators. Acta Neuropathol Commun 2023;11:79. [Crossref] [PubMed]
48. Cui JW, Li Y, Yang Y, et al. Tumor immunotherapy resistance: Revealing the mechanism of PD-1 / PD-L1-mediated tumor immune escape. Biomed Pharmacother 2024;171:116203. [Crossref] [PubMed]
49. Wang P, Chen Y, Wang C. Beyond Tumor Mutation Burden: Tumor Neoantigen Burden as a Biomarker for Immunotherapy and Other Types of Therapy. Front Oncol 2021;11:672677. [Crossref] [PubMed]
50. Gromeier M, Brown MC, Zhang G, et al. Very low mutation burden is a feature of inflamed recurrent glioblastomas responsive to cancer immunotherapy. Nat Commun 2021;12:352. [Crossref] [PubMed]
51. Sharma P, Aaroe A, Liang J, et al. Tumor microenvironment in glioblastoma: Current and emerging concepts. Neurooncol Adv 2023;5:vdad009. [Crossref] [PubMed]
52. Li K, Shi H, Zhang B, et al. Myeloid-derived suppressor cells as immunosuppressive regulators and therapeutic targets in cancer. Signal Transduct Target Ther 2021;6:362. [Crossref] [PubMed]
53. Veglia F, Sanseviero E, Gabrilovich DI. Myeloid-derived suppressor cells in the era of increasing myeloid cell diversity. Nat Rev Immunol 2021;21:485-98. [Crossref] [PubMed]
54. Yeo ECF, Brown MP, Gargett T, et al. The Role of Cytokines and Chemokines in Shaping the Immune Microenvironment of Glioblastoma: Implications for Immunotherapy. Cells 2021;10:607. [Crossref] [PubMed]
55. Chakraborty AK, Weiss A. Insights into the initiation of TCR signaling. Nat Immunol 2014;15:798-807. [Crossref] [PubMed]
56. Sridaran D, Chouhan S, Mahajan K, et al. Inhibiting ACK1-mediated phosphorylation of C-terminal Src kinase counteracts prostate cancer immune checkpoint blockade resistance. Nat Commun 2022;13:6929. [Crossref] [PubMed]
57. Chen L, Flies DB. Molecular mechanisms of T cell co-stimulation and co-inhibition. Nat Rev Immunol 2013;13:227-42. [Crossref] [PubMed]
58. Rogers BM, Smith L, Dezso Z, et al. VISTA is an activating receptor in human monocytes. J Exp Med 2021;218:e20201601. [Crossref] [PubMed]
59. Cheung PF, Yang J, Fang R, et al. Progranulin mediates immune evasion of pancreatic ductal adenocarcinoma through regulation of MHCI expression. Nat Commun 2022;13:156. [Crossref] [PubMed]
60. Shen H, Yang ES, Conry M, et al. Predictive biomarkers for immune checkpoint blockade and opportunities for combination therapies. Genes Dis 2019;6:232-46. [Crossref] [PubMed]
61. Li H, van der Merwe PA, Sivakumar S. Biomarkers of response to PD-1 pathway blockade. Br J Cancer 2022;126:1663-75. [Crossref] [PubMed]
62. Ciciola P, Cascetta P, Bianco C, et al. Combining Immune Checkpoint Inhibitors with Anti-Angiogenic Agents. J Clin Med 2020;9:675. [Crossref] [PubMed]
63. Gupta SL, Basu S, Soni V, et al. Immunotherapy: an alternative promising therapeutic approach against cancers. Mol Biol Rep 2022;49:9903-13. [Crossref] [PubMed]