A TP53-based immune prognostic model for muscle-invasive bladder cancer

Hongyan Li; Huayi Lu; Wanli Cui; Yufan Huang; Xuefei Jin

doi:10.18632/aging.202150

Research Paper Volume 13, Issue 2 pp 1929—1946

A TP53-based immune prognostic model for muscle-invasive bladder cancer

Hongyan Li^1, , Huayi Lu^2, , Wanli Cui^1, , Yufan Huang^1, , Xuefei Jin^1, ,

¹ Jilin Key Laboratory of Urologic Oncology, Department of Urology, China-Japan Union Hospital of Jilin University, Changchun, China
² The First Affiliated Hospital of USTC, Division of Life Sciences and Medicine, University of Science and Technology of China, Hefei, Anhui, China

Received: July 30, 2020 Accepted: October 10, 2020 Published: December 15, 2020

https://doi.org/10.18632/aging.202150
How to Cite

Copyright: © 2020 Li et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Muscle-invasive bladder cancer (MIBC) patients are subject to unfavorable treatment options and a high recurrence rate. The status of TP53 mutations played an essential role in the progression and the prognosis of MIBC. The present study proposed to investigate the association between TP53 mutations and immunophenotype in MIBC.

Results: We established an immune prognostic model (IPM) ground on the immune-associated genes derived from variation analysis between wild-type TP53 and mutated TP53 TCGA-MIBC patients, and validated in another cohort from GEO database. Based on IPM, we divided MIBC patients into low and high risk subgroups. The high risk MIBC patients had higher proportions of macrophages M1, and lower proportions of T cells regulatory (Tregs) and activated dendritic cells than the low risk MIBC patients. Moreover, the expression of immune checkpoints genes (PD1, CTLA4, LAG3, HAVCR2 and TIGIT) was higher in the high risk patients than the low risk patients. In clinical application, IPM exhibited better survival prediction than conventional clinical characteristics.

Conclusions: Our investigation presented practical prognostic significance for MIBC patients and displayed the overarching landscape of the immune response in the MIBC microenvironment.

Methods: Data were obtained from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO). Differentially expressed genes (DEGs) analysis between the TP53 mutated and wild-type MIBC patients was conducted. The CIBERSORT algorithm was performed to evaluate the proportion of immune cell types. Gene expression profiles from the TCGA and GEO were used as training and testing cohorts to build and validate an immune-related prognostic model (IPM). Genes in the IPM model were first screened by univariate Cox analysis, then filtered by the least absolute shrinkage and selection operator (LASSO) Cox regression. A nomogram was finally established and evaluated by combining both the IPM and other clinical factors.

Introduction

Bladder cancer is a prevalent malignant disease with 429,000 new cases and nearly 165,000 deaths worldwide annually [1, 2]. Bladder cancer develops along two different pathways: low phase non-muscle invasive bladder cancer (NMIBC) and high phase muscle-invasive bladder cancer (MIBC) [3]. Radiotherapy, combined with cisplatin-based chemotherapy, is the current standard method for high phase MIBC treatment. However, the recurrence rate of almost 70% leads MIBC patients to undergo long-term surveillance. Consequently, MIBC becomes more costly than other cancers from diagnosis to the end of life [3, 4]. Recent studies have shown that different types of tumor-infiltrating immune cells are involved in the development and prognosis of MIBC [5]. Immunotherapy has been studied as a new alternative for treating various types of cancer, especially those with an unfavorable prognosis by standard treatments [6]. Several immune-associated studies have been conducted to predict the prognosis of MIBC patients [5, 7, 8]. However, few studies have explored in detail the immune phenotype within the MIBC microenvironment and its relation with MIBC prognosis.

It has been said that no matter which orientation cancer investigation turns, TP53 gets into view. This is not only due to the essential significance in the inhibition of many human cancers but also for its leading role in various cancer-associated pathways, such as DNA repair, metabolism and antioxidant function [9]. In MIBC, patients with mutated TP53 get shorter overall survival than wild-type TP53 [10, 11]. Therefore, exploring the exact pathogenic mechanism of TP53 mutation status in MIBC and other cancers is crucial to obtain new therapeutic strategies and improve the prognosis. Remarkably, several recent studies have indicated that different immune responses are involved in the status of TP53 mutations [12, 13]. Although previous studies have found that TP53 mutations were associated with bladder cancer's clinical features, such as grade classification, cancer invasion, recurrence and poor prognosis [14]. However, their specific functions in the immune profiles of MIBC have not yet been fully elucidated. We performed a comprehensive analysis to explore the correlation between the TP53 mutation and the overall survival in the Cancer Genome Atlas (TCGA) MIBC cohort. Our results indicated that 15 immune-associated biological processes were inhibited in mutated TP53 MIBC patients. Furthermore, we identified four essential differentially expressed immune-associated genes that were posteriorly used to develop an immune prognostic model (IPM). We have shown that the proposed IPM can be employed as a useful prognostic strategy in patient management. The four identified genes can be used as potential therapeutic biomarkers for MIBC.

Results

Relationship between immunotype and TP53 somatic mutations in MIBC

As shown in Figure 1A, TP53 mutation is the universal genetic mutation in MIBC. Although many studies have explored the function of TP53 mutations and found that TP53 was associate with some of the clinical features of bladder cancer, such as grade classification, cancer invasion, recurrence and poor prognosis [14], however, the particular function on immune profiles in MIBC has not been fully elucidated. In this study, we first explored the immune-associated biological processes in MIBC regarding TP53 status using gene expression matrix, somatic mutation data of 407 MIBC samples from the TCGA cohort, and matching clinical data. The patients with (n=194) and without (n=216) TP53 mutations were used to perform gene set enrichment analysis (GSEA). These results indicated that patients from the mutated TP53 group were remarkably enriched in 506 biological processes (BH-adjusted p < 0.05) (Supplementary Table 1). In addition, 15 immune-associated biological processes were fully inhibited in this group (Figure 1B). We also explored the somatic interactions of TP53 (Figure 1C) due to cancers being caused by simultaneous mutations of multiple genes. TP53 co-occurred with FAT4, FLG, and RB1 and were mutually exclusive with FGFR3 in the TCGA MIBC cohort.

Gene set enrichment analysis between the wild-type TP53 and mutated TP53 subgroup in TCGA MIBC dataset. (A) The landscape of somatic mutations in TCGA MIBC dataset. (B) 13 inhibited immune-associated biological processes in mutated TP53 MIBC patients. (C) The landscape of somatic interactions of TP53 in TCGA MIBC dataset.

Figure 1. Gene set enrichment analysis between the wild-type TP53 and mutated TP53 subgroup in TCGA MIBC dataset. (A) The landscape of somatic mutations in TCGA MIBC dataset. (B) 13 inhibited immune-associated biological processes in mutated TP53 MIBC patients. (C) The landscape of somatic interactions of TP53 in TCGA MIBC dataset.

Identification of differentially expressed immune-associated genes in mutated and wild-type TP53 MIBC samples

To figure out the connections in TP53 mutation status with immune processes, we performed a differential expression gene analysis between wild-type and mutated TP53 groups of 1441 immune-associated genes acquired from the 15 immune-associated biological processes inhibited in the mutated TP53 group. 44 out of 1441 immune-associated genes were differentially expressed between wild-type and mutated TP53 MIBC samples (BH-adjusted p < 0.05 and absolute logFC > 1) (Supplementary Table 2).

Development of an immune prognostic model (IPM) and assessment of its predictive ability in the TCGA MIBC cohort

To assess differentially expressed immune-associated genes’ predictive ability, we performed univariate cox regression of the expression matrix of these genes. We found 9 out of 44 immune-associated genes that were notably associated with overall survival (OS) (Supplementary Table 3). To get the most significant genes associated with prognostic worth, we used cox-proportional hazards analysis based on the L1-penalized LASSO estimation. We discovered four key genes (CTSG, TREML4, KRT1 and PPBP) that occurred more than 900 times in 1000 recurrences [15, 16]. The IPM score was established according to the four key genes expression weighted by the Cox regression coefficients: IPM risk score = (0.162 × CTSG expression) + (1.167 × TREML4 expression) + (0.164 × KRT1 expression) – (0.326 × PPBP expression).

To obtain a rigid cutoff value to categorize the MIBC patients into high or low risk classes, we standardized four critical genes’ expression value with mean value = 0 and standard deviation = 1 [17]. We then used X-tile software to calculate the optimal cutoff point value (1.6) to classify patients into these two groups. In the TCGA cohort, the OS of patients in the high risk class was shorter than that of the low risk group (Figure 2A). Moreover, the high risk subgroup patients displayed a 2.30-fold higher risk (95% confidence interval (CI): 1.98–2.31, BH-adjusted p < 10^-3) than the low risk subgroup patients. The risk score and expression distribution of the four essential genes are displayed in Figure 2B. The area under the ROC curve (AUC) was used to show the predictive ability of IPM (Figure 2C).

Prognosis of IPM in TCGA MIBC cohort and meta-GEO MIBC cohort. (A, D) Kaplan-Meier survival analysis in TCGA cohort and meta-GEO cohort, OS in low risk subgroup was higher than those in the high risk subgroup. (B, E) Risk assessment in TCGA cohort and meta-GEO cohort, the correlation between the risk score (upper) and the expression of four immune-associated genes (bottom). (C, F) Time-dependent ROC curves of IPM.

Figure 2. Prognosis of IPM in TCGA MIBC cohort and meta-GEO MIBC cohort. (A, D) Kaplan-Meier survival analysis in TCGA cohort and meta-GEO cohort, OS in low risk subgroup was higher than those in the high risk subgroup. (B, E) Risk assessment in TCGA cohort and meta-GEO cohort, the correlation between the risk score (upper) and the expression of four immune-associated genes (bottom). (C, F) Time-dependent ROC curves of IPM.

Validation of the predictive ability of IPM in meta-GEO MIBC cohort

To verify the robustness of IPM constructed from the TCGA cohort, we assessed the metagene expression omnibus (meta-GEO) MIBC cohort, composed of 440 MIBC samples, with the same formula and cut off value used in the TCGA MIBC cohort. In the meta-GEO MIBC cohort, the OS of high risk group patients was lower than the low risk class (Figure 2D), as previously found in the TCGA MIBC cohort. The high risk patients displayed worse OS (HR: 2.61, 95% confidence interval (CI): 1.76–2.23, BH-adjusted p < 10^-3) than the low risk patients, indicating the applicability of the developed IPM in different populations. The risk score and expression distribution of the four essential genes are displayed in Figure 2E. The IPM obtained from AUC of the meta-GEO MIBC cohort is shown in Figure 2F. Taken together, our analysis has shown that the developed IPM has proven to be robust when faced with different datasets and, therefore, can be used for further studies.

Classification analysis of OS for IPM based on TP53 status in the TCGA MIBC cohort

As shown in Figure 3A, the phenotype of TP53 was meaningfully connected with the prognosis of patients with MIBC. To prove the relevance between the prognostic significance of IPM and TP53 status, we performed a classification analysis, separating the TCGA MIBC cohort into two groups based on TP53 status. Classification analysis demonstrated that IPM was significantly associated with OS in wild-type and mutated TP53 groups (Figure 3B, 3C). Moreover, correlation analyses showed that the risk score was significantly positively connected with OS in the wild-type TP53 group and negatively in the mutated TP53 group (Figure 3D). We performed classification analyses of several types of TP53 mutations, and we found that different TP53 mutations affect the prognosis of patients with MIBC (Figure 3E) [10, 18]. To prove the relevance between the prognostic significance of the IPM and TP53 mutations, we performed a predictive analysis of the TP53 missense mutation subgroup, which is the most common type of TP53 mutation. As expected, the TP53 missense mutation was divided into high and low risk groups by IPM (Figure 3F).

Prognosis of different TP53 mutations in TCGA MIBC cohort. (A) Kaplan-Meier survival analysis between the wild-type TP53 patients and mutated TP53 patients. (B) Kaplan-Meier survival analysis between the high and low risk subgroup in mutated TP53 MIBC patients. (C) Kaplan-Meier survival analysis between the high and low risk subgroup in wild-type TP53 MIBC patients. (D) Correlation analysis between risk score and survival time according to TP53 status. (E) Kaplan-Meier survival analysis among the different types of TP53mutations. (F) Kaplan-Meier survival analysis between the high and low risk subgroup in TP53 MIBC missense mutation patients.

Figure 3. Prognosis of different TP53 mutations in TCGA MIBC cohort. (A) Kaplan-Meier survival analysis between the wild-type TP53 patients and mutated TP53 patients. (B) Kaplan-Meier survival analysis between the high and low risk subgroup in mutated TP53 MIBC patients. (C) Kaplan-Meier survival analysis between the high and low risk subgroup in wild-type TP53 MIBC patients. (D) Correlation analysis between risk score and survival time according to TP53 status. (E) Kaplan-Meier survival analysis among the different types of TP53mutations. (F) Kaplan-Meier survival analysis between the high and low risk subgroup in TP53 MIBC missense mutation patients.

Immune landscape comparison between the low and high risk MIBC patients

To assess the variation of 22 immune infiltrated cell types between low and high risk MIBC patients, we performed an immune infiltration analysis using CIBERSORT analytical tool with the LM22 signature matrix [19]. The immune landscape of the 407 TCGA MIBC cohort was summarized by representing the proportions of the 22 immune cell types within and between low and high risk groups (Figure 4A). The proportions of diverse subtypes of tumor-infiltrating immune cells were slightly moderately correlated (Figure 4B). The high risk MIBC patients had significantly higher proportions of macrophages M1, while expressively lower proportions of T cells regulatory (Tregs) and activated dendritic cells than the low risk MIBC patients (BH-adjusted p < 0.05) (Figure 4C). Hence, our results indicated that aberrant immune cell infiltration and the heterogeneity of tumor infiltration in MIBC might serve as prognostic indicators and targets for immunotherapy and have significant clinical implications.

Comparison of immune infiltration landscape between the high and low risk MIBC patients. (A) Immune cell proportions between the high and low risk MIBC patients. (B) Correlation matrix of 22 types of immune cell proportions. (C) Significantly difference of immune cells between the high and low risk MIBC patients.

Figure 4. Comparison of immune infiltration landscape between the high and low risk MIBC patients. (A) Immune cell proportions between the high and low risk MIBC patients. (B) Correlation matrix of 22 types of immune cell proportions. (C) Significantly difference of immune cells between the high and low risk MIBC patients.

Immune checkpoint inhibition has been a way explored by many novel anti-tumor agents, which have been used successfully in different types of cancer [20]. In bladder cancer, drugs targeting immune checkpoints genes have also been indicated to play anti-cancer function by reversing tumor immunosuppressive effects [21]. Therefore, we estimated the correlation between the risk score of MIBC patients and expression level of vital immune checkpoints genes (PD1, CTLA4, LAG3, HAVCR2 and TIGIT) (Supplementary Table 4). We demonstrated that the risk score of MIBC patients was significantly associated with the expression of PD1, CTLA4, LAG3, HAVCR2 and TIGIT (BH-adjusted p < 0.05) (Figure 5A). Furthermore, we showed that the expression level of PD1, CTLA4, LAG3, HAVCR2 and TIGIT genes in the high risk MIBC subgroup was significantly higher than in the low risk subgroup (BH-adjusted p < 0.01) (Figure 5B–5F). These data suggest that the bleak prognosis of high risk MIBC patients may be partially due to the immunosuppressive microenvironment.

Enrichment analysis of IPM. (A) Correlation between the risk score and the expression level of immune checkpoints genes. (B–F) The expression level of immune checkpoints genes in the high and low risk MIBC patients. (G) Heatmap plot of immune-associated genes that were differentially expressed between the high and low risk MIBC patients. (H) Enrichment of biological processes for immune-associated genes that are shown in the circular plot. (I) Enrichment of KEGG pathways for immune-associated genes that are shown in the sankey plot.

Figure 5. Enrichment analysis of IPM. (A) Correlation between the risk score and the expression level of immune checkpoints genes. (B–F) The expression level of immune checkpoints genes in the high and low risk MIBC patients. (G) Heatmap plot of immune-associated genes that were differentially expressed between the high and low risk MIBC patients. (H) Enrichment of biological processes for immune-associated genes that are shown in the circular plot. (I) Enrichment of KEGG pathways for immune-associated genes that are shown in the sankey plot.

Variation of biological processes and pathways in low and high risk MIBC subgroup patients

We filtered the 44 immune-associated genes differentially expressed between low and high risk subgroups of MIBC patients by risk score association using the correlation analysis method. This analysis identified 25 immune-associated genes as risk score associated (Pearson correlation coefficient >0.2 and BH-adjusted p < 0.05; Figure 5G). GO and KEGG analyses were used to establish the inherent biological functions and pathways involved in these genes (FDR < 0.0001 and FDR < 0.001, respectively; Figures 5H, 5I) (Supplementary Tables 5 and 6). The 25 immune-associated genes related to the risk score in the TCGA MIBC cohort were principally enriched in the humoral and innate immune responses and T cell receptor signaling pathway (Figures 5H, 5I).

IPM exhibits a better survival prediction than conventional clinical characteristics

To compare the prognostic value of IPM with conventional clinical characteristics in the TCGA MIBC cohort, we performed univariate and multivariate Cox regression analyses involving five conventional clinical characteristics (age, gender, weight, subtype, and pathological stage). The univariate Cox regression analysis showed that IPM was independent of the prognostic aspect, accordingly proving its robustness to predict MIBC prognosis independently (Figure 6A). Furthermore, the multivariate Cox regression analysis demonstrated that IPM was meaningly associated with the survival time (BH-adjusted p < 0.01) and the highest median risk score (HR = 2.22, 95% CI = 1.38 – 3.56). In addition to the Cox regression analysis, we also estimated the c-index between IPM and conventional clinical characteristics. The c-index value of IPM (0.7181) was higher than other conventional clinical characteristics (0.1821 – 0.6438) (Figure 6B). Taken together, these results demonstrated that IPM has a better survival prediction than conventional clinical characteristics.

The connection between IPM and conventional clinical characteristics. (A) Univariate and multivariate regression analysis of IPM and clinical characteristics in prognostic value. Blue displays no statistical significance, and red displays statistical significance. (B) Comparison of c-index among different clinical characteristics. (C) Nomogram for predicting the probability of 1, 3, and 5 years OS for MIBC patients. (D) Calibration plot of the nomogram for predicting the probability of OS at 1, 3, and 5 years (red lines); dash lines indicate the actual probability. The vertical mark lines in the top x-axis (marginal rugs) stand for the distribution of samples in each fitting curve.

Figure 6. The connection between IPM and conventional clinical characteristics. (A) Univariate and multivariate regression analysis of IPM and clinical characteristics in prognostic value. Blue displays no statistical significance, and red displays statistical significance. (B) Comparison of c-index among different clinical characteristics. (C) Nomogram for predicting the probability of 1, 3, and 5 years OS for MIBC patients. (D) Calibration plot of the nomogram for predicting the probability of OS at 1, 3, and 5 years (red lines); dash lines indicate the actual probability. The vertical mark lines in the top x-axis (marginal rugs) stand for the distribution of samples in each fitting curve.

Construction and validation of the individualized nomogram on IPM

To predict the prognosis of MIBC patients for clinicians, we used independently predictors (age, subtype, pathologic stage and IPM) that were identified in the multivariate Cox regression analysis (BH-adjusted p < 0.05) to construct the nomogram. As shown in Figure 6C, the IPM presented more risk points than other predictors. The c-index of the nomogram was 0.72 (95%CI, 0.61 to 0.86) using 1000 bootstrap replicates. The calibration curve displayed good agreement between the prediction and the investigation (Figure 6D). Except for the probability of 1-year survival, others were significantly lower than that of the nomogram (BH-adjusted p

Discussion

TP53 mutations are related to high grade, invasive tumor, low recurrence and adverse clinical outcomes in bladder cancer [14]. Recent studies have found that TP53 mutation can increase the gene expression involved in immune checkpoints, activate T-effector cells, and elevate interferon-γ in lung adenocarcinoma [12]. Moreover, co-occurring TP53/KRAS mutation displayed more benefit from PD-1 inhibitors [12]. TP53 also was reported to be a predictor in immunotherapy of PD-1 in lung carcinomatosis [22]. Regarding bladder cancer, it has been demonstrated that anti-PD-1 antibodies enhance radiotherapy-induced anti-tumor immunity. However, the molecular mechanism of TP53 mutations in MIBC immunophenotype regulation and MIBC's prognosis is unknown. Therefore, the immune-associated effects of TP53 mutations in bladder cancer must be studied.

Furthermore, the development of advanced immune-associated prognostic models would be useful to identify biomarkers, evaluate the immune state of MIBC patients, and classify them to enhance the effectiveness of immunotherapy. In recent years, tumor immune signatures have been established and identified in various cancers [23–25]. Even though some studies have attempted to illustrate immune signatures in bladder cancer [26, 27], the essence of the local immune landscape in MIBC prognosis and prediction has not been fully elucidated. In the present study, we explored the character of TP53 mutations in the modulation of immune signature in MIBC. The GSEA analysis showed that the mutated TP53 group had a significantly lower immune phenotype than wild-type TP53. Moreover, differential expression analysis of these immune-associated genes showed that 15 immune-associated biological processes were all inhibited in the mutated TP53 group. Subsequently, a cox-proportional hazards analysis based on the L1-penalized LASSO estimation pointed four key genes (CTSG, TREML4, KRT1 and PPBP) used to construct a novel immune prognostic model (IPM) to predict the prognosis of MIBC patients.

These four key genes are involved in various cancer-related immune processes. Cathepsin G (CTSG), an azurophil granule protease, can lead to breast cancer cell migration [28, 29] and can enter tumor endosomes by binding to a cell surface receptor [30]. Previous studies have shown that CTSG degraded MHC I on the human glioblastoma cell surface of primary immune cells. CTSG activity has been presented as a new way to glioblastoma treatment [31]. Besides, high CTSG expression was correlated with poor outcomes in treating acute myeloid leukemia (AML) [32]. The triggering receptor expressed on myeloid cells 4 (TREML4) is a Ig superfamily member and rarely reported in cancer [33]. However, we believe that TREML4 could be a new immunotherapeutic target because a TREML4 deficiency causes a fail in INF-I production by macrophages due to a decrease of phosphorylation level in the signal transducer and activator of transcription 1 (STAT1). Moreover, previous studies have found that TREM1 dramatically promotes proliferation and invasion in Hepatocellular Carcinoma (HCC) and TREM1 expression has been associated with poor survival in HCC patients [34]. In addition, high TREM2 expression was inversely associated with unfavorable prognosis in gastric cancer and, therefore, could be a useful prognostic biomarker [35]. Keratin 1 (KRT1) is a member of the keratin family, which is involved in protein binding [36], carbohydrate-binding [37] and structural constituent of the epidermis [38]. Recently, it was found that KRT1 was highly expressed in breast cancer cells, and, thus, pointed as a new marker for breast cancer targeting [39]. In addition, KRT1 was significantly overexpressed in colon cancer, especially in its later phase [39]. Finally, pro-platelet basic protein (PPBP), also known as CXCL7, is a biomarker of several cancer types, such as renal cell cancer [40], lung cancer [41] and colorectal cancer [42]. PPBP has been reported to be up-regulated in the peripheral blood of early-stage renal cell carcinoma patients [40]. PPBP accelerated the development of renal cell carcinoma by promoting cell proliferation [43]. In colorectal cancer, the expression of PPBP was significantly correlated with sex, TNM and T stages [42].

We demonstrated by different approaches that our IPM exhibited better survival prediction than traditional clinical features. Subsequently, we performed a comprehensive evaluation that integrated the IPM with conventional clinical characteristics (age, gender, weight, subtype, and pathological stage). The calibration curve displayed good agreement between the prediction and the clinical characteristics evaluated for 1, 3 and 5 years OS. Our main superiority is that it presents a complementary perspective on individual tumors and establishes a unique scoring method for MIBC patients. Therefore, our nomogram could be used as a promising tool for clinicians in the future.

The stage of elimination is an upgraded model of cancer immunosurveillance. The innate and adaptive immune systems co-operate to find the existence of the progression of tumors and devastate it before it becomes worse [44]. Like Type I IFNs, some danger signals can induce the growth of tumors, motivate dendritic cells and promote the adaptive immune responses [45]. In the equilibrium stage, the adaptive immune system restrains tumor cell outgrowth and forms the tumor cells’ immunogenicity. IL-12, IFN-γ, CD4⁺ and CD8⁺ T cells are responsible for retaining the occult tumor cells [44, 46].

The cancer immunoediting process involves three stages: elimination, equilibrium and escape [44]. In the escape stage, tumor cells that have passed through the former two steps and have obtained the ability to escape from immune recognition become visible and progressively growing tumors. Therefore, the occurrence of tumor escape is based on the establishment of an immunosuppressive status within the tumor microenvironment [44, 47]. Regulatory T cells (Tregs) are one of the main types of immunosuppressive cells responsible for restraining host-protective anti-tumor responses. Activated Tregs can express PD1, PDL1 and CTLA4 to inhibit tumor-specific T lymphocytes functions [44, 48]. We compared the immune infiltration of 22 immune cell types between low and high risk MIBC patients to explore the immune mechanisms and evaluate the reach of the proposed IPM as cancer immunotherapy. The results showed that the high risk MIBC patients had more macrophages. The patients with (n=194) and without (n=216) TP53 mutations, and fewer Tregs and activated dendritic cells than the low risk MIBC patients. Furthermore, the gene expression levels of PD1, CTLA4, LAG3, HAVCR2 and TIGIT in the high risk MIBC subgroup were significantly higher than those in the low risk subgroup. Therefore, the risk score obtained from IPM was consistent with the ability of immune infiltration to decide the expressed value of immune checkpoint genes. This suggests that the inferior prognosis of the high risk patients may have resulted from a more excellent immunosuppressive environment and an increased expression level of immune checkpoint genes. Consequently, our results also suggest that immune checkpoint gene inhibitors should be potentially more effective in high risk MIBC patients, which would result in a considerable improvement in prognosis.

Conclusion

The appointments presented here can present some limitations due to being based on retrospective data and, therefore, they should be further validated by prospective studies. Furthermore, the four key immune-associated genes used to construct the IPM should also be detected in experimental studies to ensure their clinical application. On the other hand, the present work provides new insights into the MIBC immune microenvironment and immune-associated therapies. For the first time, it was proposed an IPM based on TP53 mutations and four immune-associated genes. The proposed IPM presented a significantly effective prognostic for MIBC patients and illustrated an overarching landscape of immune response in the MIBC microenvironment. Remarkably, the development and validation of the IPM presented an immunological perspective to elucidate the mechanisms in the clinical outcomes of MIBC and potentially could be used as a reference for the study of other types of cancer.

Materials and Methods

Cancer Genome Atlas (TCGA) data acquisition and processing

The gene expression matrix, somatic mutation data of 407 MIBC samples and their matching clinical features were acquired from the Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/repository). RNA-seq count data were downloaded using Illumina HiSeq platforms and annotated as Homo_sapiens. GRCh38.91.chr.gtf file (http://asia.ensembl.org/index.html). The function of variances stabilizing transformation (VST) normalization method of the DESeq2 package in R software was used to normalize the expression data [49]. The mean gene expression level was employed when multiple gene symbols were located. Moreover, we removed the genes that sum gene expression value to be less than 100 to filter out the low abundance data.

Gene Expression Omnibus (GEO) data acquisition and processing

The four gene expression matrix of microarray data from GSE13507 based on platform GPL6102 (257 samples used in this study), GSE32549 based on platform GPL6947 (132 samples used in this study), GSE48075 based on platform GPL6947 (143 samples used in this study), and GSE48276 based on platform GPL14951 (117 samples used in this study) were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The mean gene expression level was employed when multiple gene symbols were located. Genes were removed when its sum gene expression value was less than 100 to filter out the low abundance data. Moreover, the clinical data of four datasets with survival details were combined into the meta-GEO MIBC cohort (n=440) from 649 GEO patients to confirm the immune prognostic model. Combat method in sva R package was applied to remove the batch effects [50, 51].

Gene Set Enrichment Analysis (GSEA)

We divided the TCGA cohort into two groups, wild-type TP53 (n = 216) and mutated TP53 (n = 194), to show the variations in immunological pathways and related genes using GSEA methods [52]. The MSigDB gene sets file (c5.bp.v7.0.symbols.gmt) was chosen as the reference gene set with permutations of 10⁴. FDR < 0.05 was appointed as the threshold.

Differentially expressed genes (DEGs) and functional analysis

The R package DESeq2 was used to filter out the differentially expressed genes (DEGs) between the wild-type and mutated TP53 groups. Absolute log₂(fold change) > 2 and adjusted p values < 0.01 were set as the cut-off criteria to indicate significant statistically difference [49]. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analyses were performed using the clusterProfiler R package to explore the significant biological processes, highlight the pathways associated with DEGs, and assess the biological implications of the prognostic model [53]. The notable biological processes and pathways were presented using GOplot [54] and ggalluvial (version 0.10.0. https://CRAN.R-project.org/package=ggalluvial) R packages, respectively.

Immune prognostic model (IPM) development and validation

To predict the prognosis of patients with MIBC, we constructed an immune prognostic model (IPM) as previously described [55]. We integrated 407 MIBC samples expression profiles, mutation data and survival information from TCGA MIBC samples to subsequent analyses. The expression matrix of the DEGs in immune-associated genes sets of GSEA from 407 MIBC was analyzed using univariate cox regression analysis. The significant genes with BH-adjusted p value < 0.01 were subject to the least absolute shrinkage and selection operator (LASSO) analysis with L1-penalty. LASSO is a widely used method for dealing with the very high dimensional space of predictors such as gene expression profiles [56]. Thus, the critical immune-associated genes, significant in univariate cox regression analysis, were screened out by the LASSO method. Finally, a comparatively small portion of non-zero weight genes remained. Most of the potential indicators were reduced to zero. Hence, we decreased the number of immune genes using LASSO-penalized cox regression. In this study, we performed LASSO-cox analysis using glmnet R package, and picked up immune-associated genes that appeared more than 900 times in 1000 repetitions [15, 57]. X-tile 3.6.1 software (Yale University, New Haven, CT, USA) was employed to define the best cutoff for categorizing low or high risk MIBC patients. The predictive ability of the IPM was assessed by the log-rank test and Kaplan-Meier survival analysis.

Estimation of immune cell type fractions

CIBERSORT analytical tool [19] was employed to quantify the immune cell distributions in wild-type TP53 and mutated TP53 MIBC. The LM22 gene signature was designed to evaluate the possibility of leukocyte deconvolution from bulk tumors. LM22 signature matrix file contains 547 genes and enables highly sensitive and particular distinction of 22 human hematopoietic cell phenotypes, including T cell subtypes, B cell subtypes, NK cell subtypes, plasma cells subtypes and myeloid subtypes. For every sample, the sum of all estimates of 22 immune cell subtype fractions was equivalent to 1.

Comparison between IPM and traditional clinical characteristics

We performed univariate and multivariate cox regression analyses for MIBC patients to compare the predictive ability of IPM and traditional clinical characteristics. The traditional clinical characteristics used here were obtained from survival information of 407 MIBC patients, including age, gender, weight, pathological stage and diagnosis type.

Assessment of nomogram performance

The nomogram was used to show the predicted probability of survival for 1 year, 3 and 5 years, based on the outcome of multivariate analysis. The calibration curve and concordance index (c-index), which were generated by the rms (Version: 5.1–3.1; https://CRAN.R-project.org/package=rms) R package [58], decided the predictive exactness and discriminative capacity of a nomogram. The calibration curve indicated the distinction between the actual overall survival rate and the nomogram's predicted probability. A calibration curve closer to the diagonal dotted line suggests a better predictive effect. The c-index is mainly used to assess the predictive power of the model. The model impact can be equivalent to the area under the ROC curve (AUC). We calculated c-index by a bootstrap approach with resampling 1000 times [59]. For the statistical tests analyzed in this study, BH-adjusted p < 0.05 was used as a threshold to indicate statistically significant differences.

Supplementary Materials

Supplementary Table 1

Supplementary Tables 2, 3 and 4

Supplementary Table 5

Supplementary Table 6

Abbreviations

MIBC: Muscle-invasive bladder cancer; IPM: immune prognostic model; Tregs: T cells regulatory; NMIBC: non-muscle-invasive bladder cancer; TCGA: the Cancer Genome Atlas; VST: variance stabilizing transformation; GEO: Gene Expression Omnibus; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; LASSO: least absolute shrinkage and selection operator; AUC: area under the ROC curve; GSEA: gene set enrichment analysis; CTSG: Cathepsin G; AML: acute myeloid leukemia; TREML4: Triggering receptor expressed on myeloid cells 4; STAT1: signal transducer and activator of transcription 1.

Author Contributions

XJ conceived and designed the study. HL collected and analyzed the data. HL and HYL wrote the manuscript. WC and YH revised the manuscript. All authors read and approved the manuscript.

Conflicts of Interest

These authors declare no conflicts of interest.

References

1. Dyrskjøt L, Zieger K, Real FX, Malats N, Carrato A, Hurst C, Kotwal S, Knowles M, Malmström PU, de la Torre M, Wester K, Allory Y, Vordos D, et al. Gene expression signatures predict outcome in non-muscle-invasive bladder carcinoma: a multicenter validation study. Clin Cancer Res. 2007; 13:3545–51. https://doi.org/10.1158/1078-0432.CCR-06-2940 [PubMed]
2. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, Parkin DM, Forman D, Bray F. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer. 2015; 136:E359–86. https://doi.org/10.1002/ijc.29210 [PubMed]
3. Botteman MF, Pashos CL, Redaelli A, Laskin B, Hauser R. The health economics of bladder cancer: a comprehensive review of the published literature. Pharmacoeconomics. 2003; 21:1315–30. https://doi.org/10.1007/BF03262330 [PubMed]
4. Kaufman DS, Shipley WU, Feldman AS. Bladder cancer. Lancet. 2009; 374:239–49. https://doi.org/10.1016/S0140-6736(09)60491-8 [PubMed]
5. Liu Z, Zhu Y, Xu L, Zhang J, Xie H, Fu H, Zhou Q, Chang Y, Dai B, Xu J. Tumor stroma-infiltrating mast cells predict prognosis and adjuvant chemotherapeutic benefits in patients with muscle invasive bladder cancer. Oncoimmunology. 2018; 7:e1474317. https://doi.org/10.1080/2162402X.2018.1474317 [PubMed]
6. Long J, Lin J, Wang A, Wu L, Zheng Y, Yang X, Wan X, Xu H, Chen S, Zhao H. PD-1/PD-L blockade in gastrointestinal cancers: lessons learned and the road toward precision immunotherapy. J Hematol Oncol. 2017; 10:146. https://doi.org/10.1186/s13045-017-0511-2 [PubMed]
7. Efstathiou JA, Mouw KW, Gibb EA, Liu Y, Wu CL, Drumm MR, da Costa JB, du Plessis M, Wang NQ, Davicioni E, Feng FY, Seiler R, Black PC, et al. Impact of immune and stromal infiltration on outcomes following bladder-sparing trimodality therapy for muscle-invasive bladder cancer. Eur Urol. 2019; 76:59–68. https://doi.org/10.1016/j.eururo.2019.01.011 [PubMed]
8. Song BN, Kim SK, Mun JY, Choi YD, Leem SH, Chu IS. Identification of an immunotherapy-responsive molecular subtype of bladder cancer. EBioMedicine. 2019; 50:238–45. https://doi.org/10.1016/j.ebiom.2019.10.058 [PubMed]
9. Levine AJ, Oren M. The first 30 years of p53: growing ever more complex. Nat Rev Cancer. 2009; 9:749–58. https://doi.org/10.1038/nrc2723 [PubMed]
10. Nassar AH, Umeton R, Kim J, Lundgren K, Harshman L, Van Allen EM, Preston M, Dong F, Bellmunt J, Mouw KW, Choueiri TK, Sonpavde G, Kwiatkowski DJ. Mutational analysis of 472 urothelial carcinoma across grades and anatomic sites. Clin Cancer Res. 2019; 25:2458–70. https://doi.org/10.1158/1078-0432.CCR-18-3147 [PubMed]
11. Dadhania V, Zhang M, Zhang L, Bondaruk J, Majewski T, Siefker-Radtke A, Guo CC, Dinney C, Cogdell DE, Zhang S, Lee S, Lee JG, Weinstein JN, et al. Meta-analysis of the luminal and basal subtypes of bladder cancer and the identification of signature immunohistochemical markers for clinical use. EBioMedicine. 2016; 12:105–17. https://doi.org/10.1016/j.ebiom.2016.08.036 [PubMed]
12. Dong ZY, Zhong WZ, Zhang XC, Su J, Xie Z, Liu SY, Tu HY, Chen HJ, Sun YL, Zhou Q, Yang JJ, Yang XN, Lin JX, et al. Potential predictive value of TP53 and KRAS mutation status for response to PD-1 blockade immunotherapy in lung adenocarcinoma. Clin Cancer Res. 2017; 23:3012–24. https://doi.org/10.1158/1078-0432.CCR-16-2554 [PubMed]
13. Smal MP, Rolevich AI, Polyakov SL, Krasny SA, Goncharova RI. FGFR3 and TP53 mutations in a prospective cohort of Belarusian bladder cancer patients. Exp Oncol. 2014; 36:246–51. [PubMed]
14. Esrig D, Elmajian D, Groshen S, Freeman JA, Stein JP, Chen SC, Nichols PW, Skinner DG, Jones PA, Cote RJ. Accumulation of nuclear p53 and tumor progression in bladder cancer. N Engl J Med. 1994; 331:1259–64. https://doi.org/10.1056/NEJM199411103311903 [PubMed]
15. Xu RH, Wei W, Krawczyk M, Wang W, Luo H, Flagg K, Yi S, Shi W, Quan Q, Li K, Zheng L, Zhang H, Caughey BA, et al. Circulating tumour DNA methylation markers for diagnosis and prognosis of hepatocellular carcinoma. Nat Mater. 2017; 16:1155–61. https://doi.org/10.1038/nmat4997 [PubMed]
16. Tibshirani R. Regression Shrinkage and Selection Via the Lasso. Journal of the Royal Statistical Society. 2011; 58:267–288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x
17. Shen S, Wang G, Zhang R, Zhao Y, Yu H, Wei Y, Chen F. Development and validation of an immune gene-set based prognostic signature in ovarian cancer. EBioMedicine. 2019; 40:318–26. https://doi.org/10.1016/j.ebiom.2018.12.054 [PubMed]
18. Pietzak EJ, Bagrodia A, Cha EK, Drill EN, Iyer G, Isharwal S, Ostrovnaya I, Baez P, Li Q, Berger MF, Zehir A, Schultz N, Rosenberg JE, et al. Next-generation sequencing of nonmuscle invasive bladder cancer reveals potential biomarkers and rational therapeutic targets. Eur Urol. 2017; 72:952–59. https://doi.org/10.1016/j.eururo.2017.05.032 [PubMed]
19. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–57. https://doi.org/10.1038/nmeth.3337 [PubMed]
20. Walshaw RC, Honeychurch J, Illidge TM, Choudhury A. The anti-PD-1 era - an opportunity to enhance radiotherapy for patients with bladder cancer. Nat Rev Urol. 2018; 15:251–59. https://doi.org/10.1038/nrurol.2017.172 [PubMed]
21. Morales A, Eidinger D, Bruce AW. Intracavitary bacillus calmette-guerin in the treatment of superficial bladder tumors. J Urol. 2017; 197:S142–45. https://doi.org/10.1016/j.juro.2016.10.101 [PubMed]
22. Skoulidis F, Goldberg ME, Greenawalt DM, Hellmann MD, Awad MM, Gainor JF, Schrock AB, Hartmaier RJ, Trabucco SE, Gay L, Ali SM, Elvin JA, Singal G, et al. STK11/LKB1 mutations and PD-1 inhibitor resistance in KRAS-mutant lung adenocarcinoma. Cancer Discov. 2018; 8:822–35. https://doi.org/10.1158/2159-8290.CD-18-0099 [PubMed]
23. Jiang Y, Zhang Q, Hu Y, Li T, Yu J, Zhao L, Ye G, Deng H, Mou T, Cai S, Zhou Z, Liu H, Chen G, et al. ImmunoScore signature: a prognostic and predictive tool in gastric cancer. Ann Surg. 2018; 267:504–13. https://doi.org/10.1097/SLA.0000000000002116 [PubMed]
24. Li Y, Lu Z, Che Y, Wang J, Sun S, Huang J, Mao S, Lei Y, Chen Z, He J. Immune signature profiling identified predictive and prognostic factors for esophageal squamous cell carcinoma. Oncoimmunology. 2017; 6:e1356147. https://doi.org/10.1080/2162402X.2017.1356147 [PubMed]
25. Lu X, Jiang L, Zhang L, Zhu Y, Hu W, Wang J, Ruan X, Xu Z, Meng X, Gao J, Su X, Yan F. Immune signature-based subtypes of cervical squamous cell carcinoma tightly associated with human papillomavirus type 16 expression, molecular features, and clinical outcome. Neoplasia. 2019; 21:591–601. https://doi.org/10.1016/j.neo.2019.04.003 [PubMed]
26. Pfannstiel C, Strissel PL, Chiappinelli KB, Sikic D, Wach S, Wirtz RM, Wullweber A, Taubert H, Breyer J, Otto W, Worst T, Burger M, Wullich B, et al, and BRIDGE Consortium, Germany. The tumor immune microenvironment drives a prognostic relevance that correlates with bladder cancer subtypes. Cancer Immunol Res. 2019; 7:923–38. https://doi.org/10.1158/2326-6066.CIR-18-0758 [PubMed]
27. Pichler R, Fritz J, Zavadil C, Schäfer G, Culig Z, Brunner A. Tumor-infiltrating immune cell subpopulations influence the oncologic outcome after intravesical bacillus calmette-guérin therapy in bladder cancer. Oncotarget. 2016; 7:39916–30. https://doi.org/10.18632/oncotarget.9537 [PubMed]
28. Morimoto-Kamata R, Yui S. Insulin-like growth factor-1 signaling is responsible for cathepsin G-induced aggregation of breast cancer MCF-7 cells. Cancer Sci. 2017; 108:1574–83. https://doi.org/10.1111/cas.13286 [PubMed]
29. Morimoto-Kamata R, Mizoguchi S, Ichisugi T, Yui S. Cathepsin G induces cell aggregation of human breast cancer MCF-7 cells via a 2-step mechanism: catalytic site-independent binding to the cell surface and enzymatic activity-dependent induction of the cell aggregation. Mediators Inflamm. 2012; 2012:456462. https://doi.org/10.1155/2012/456462 [PubMed]
30. Gregory AD, Hale P, Perlmutter DH, Houghton AM. Clathrin pit-mediated endocytosis of neutrophil elastase and cathepsin G by cancer cells. J Biol Chem. 2012; 287:35341–50. https://doi.org/10.1074/jbc.M112.385617 [PubMed]
31. Palesch D, Wagner J, Meid A, Molenda N, Sienczyk M, Burkhardt J, Münch J, Prokop L, Stevanovic S, Westhoff MA, Halatsch ME, Wirtz CR, Zimecki M, Burster T. Cathepsin G-mediated proteolytic degradation of MHC class I molecules to facilitate immune detection of human glioblastoma cells. Cancer Immunol Immunother. 2016; 65:283–91. https://doi.org/10.1007/s00262-016-1798-5 [PubMed]
32. Alatrash G, Garber HR, Zhang M, Sukhumalchandra P, Qiu Y, Jakher H, Perakis AA, Becker L, Yoo SY, Dwyer KC, Coombes K, Talukder AH, John LS, et al. Cathepsin G is broadly expressed in acute myeloid leukemia and is an effective immunotherapeutic target. Leukemia. 2017; 31:234–37. https://doi.org/10.1038/leu.2016.249 [PubMed]
33. Sen SK, Boelte KC, Barb JJ, Joehanes R, Zhao X, Cheng Q, Adams L, Teer JK, Accame DS, Chowdhury S, Singh LN, Kavousi M, Peyser PA, et al, and NISC Comparative Sequencing Program, and CHARGE Consortium. Integrative DNA, RNA, and protein evidence connects TREML4 to coronary artery calcification. Am J Hum Genet. 2014; 95:66–76. https://doi.org/10.1016/j.ajhg.2014.06.003 [PubMed]
34. Duan M, Wang ZC, Wang XY, Shi JY, Yang LX, Ding ZB, Gao Q, Zhou J, Fan J. TREM-1, an inflammatory modulator, is expressed in hepatocellular carcinoma cells and significantly promotes tumor progression. Ann Surg Oncol. 2015; 22:3121–29. https://doi.org/10.1245/s10434-014-4191-7 [PubMed]
35. Zhang X, Wang W, Li P, Wang X, Ni K. High TREM2 expression correlates with poor prognosis in gastric cancer. Hum Pathol. 2018; 72:91–99. https://doi.org/10.1016/j.humpath.2017.10.026 [PubMed]
36. Mahdi F, Shariat-Madar Z, Todd RF 3rd, Figueroa CD, Schmaier AH. Expression and colocalization of cytokeratin 1 and urokinase plasminogen activator receptor on endothelial cells. Blood. 2001; 97:2342–50. https://doi.org/10.1182/blood.v97.8.2342 [PubMed]
37. Steinert PM, Marekov LN. The proteins elafin, filaggrin, keratin intermediate filaments, loricrin, and small proline-rich proteins 1 and 2 are isodipeptide cross-linked components of the human epidermal cornified cell envelope. J Biol Chem. 1995; 270:17702–11. https://doi.org/10.1074/jbc.270.30.17702 [PubMed]
38. Collard CD, Montalto MC, Reenstra WR, Buras JA, Stahl GL. Endothelial oxidative stress activates the lectin complement pathway: role of cytokeratin 1. Am J Pathol. 2001; 159:1045–54. https://doi.org/10.1016/S0002-9440(10)61779-8 [PubMed]
39. Soudy R, Etayash H, Bahadorani K, Lavasanifar A, Kaur K. Breast cancer targeting peptide binds keratin 1: a new molecular marker for targeted drug delivery to breast cancer. Mol Pharm. 2017; 14:593–604. https://doi.org/10.1021/acs.molpharmaceut.6b00652 [PubMed]
40. Kinouchi T, Uemura M, Wang C, Ishizuya Y, Yamamoto Y, Hayashi T, Matsuzaki K, Nakata W, Yoshida T, Jingushi K, Kawashima A, Ujike T, Nagahara A, et al. Expression level of CXCL7 in peripheral blood cells is a potential biomarker for the diagnosis of renal cell carcinoma. Cancer Sci. 2017; 108:2495–502. https://doi.org/10.1111/cas.13414 [PubMed]
41. Du Q, Li E, Liu Y, Xie W, Huang C, Song J, Zhang W, Zheng Y, Wang H, Wang Q. CTAPIII/CXCL7: a novel biomarker for early diagnosis of lung cancer. Cancer Med. 2018; 7:325–35. https://doi.org/10.1002/cam4.1292 [PubMed]
42. Li L, Zhang L, Tian Y, Zhang T, Duan G, Liu Y, Yin Y, Hua D, Qi X, Mao Y. Serum chemokine CXCL7 as a diagnostic biomarker for colorectal cancer. Front Oncol. 2019; 9:921. https://doi.org/10.3389/fonc.2019.00921 [PubMed]
43. Grépin R, Guyot M, Giuliano S, Boncompagni M, Ambrosetti D, Chamorey E, Scoazec JY, Negrier S, Simonnet H, Pagès G. The CXCL7/CXCR1/2 axis is a key driver in the growth of clear cell renal cell carcinoma. Cancer Res. 2014; 74:873–83. https://doi.org/10.1158/0008-5472.CAN-13-1267 [PubMed]
44. Schreiber RD, Old LJ, Smyth MJ. Cancer immunoediting: integrating immunity’s roles in cancer suppression and promotion. Science. 2011; 331:1565–70. https://doi.org/10.1126/science.1203486 [PubMed]
45. Matzinger P. Tolerance, danger, and the extended family. Annu Rev Immunol. 1994; 12:991–1045. https://doi.org/10.1146/annurev.iy.12.040194.005015 [PubMed]
46. Loeser S, Loser K, Bijker MS, Rangachari M, van der Burg SH, Wada T, Beissert S, Melief CJ, Penninger JM. Spontaneous tumor rejection by cbl-b-deficient CD8+ T cells. J Exp Med. 2007; 204:879–91. https://doi.org/10.1084/jem.20061699 [PubMed]
47. Radoja S, Rao TD, Hillman D, Frey AB. Mice bearing late-stage tumors have normal functional systemic T cell responses in vitro and in vivo. J Immunol. 2000; 164:2619–28. https://doi.org/10.4049/jimmunol.164.5.2619 [PubMed]
48. Dunn GP, Bruce AT, Ikeda H, Old LJ, Schreiber RD. Cancer immunoediting: from immunosurveillance to tumor escape. Nat Immunol. 2002; 3:991–98. https://doi.org/10.1038/ni1102-991 [PubMed]
49. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010; 11:R106. https://doi.org/10.1186/gb-2010-11-10-r106 [PubMed]
50. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical bayes methods. Biostatistics. 2007; 8:118–27. https://doi.org/10.1093/biostatistics/kxj037 [PubMed]
51. Leek JT, Storey JD. Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLoS Genet. 2007; 3:1724–35. https://doi.org/10.1371/journal.pgen.0030161 [PubMed]
52. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102 [PubMed]
53. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
54. Walter W, Sánchez-Cabo F, Ricote M. GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics. 2015; 31:2912–14. https://doi.org/10.1093/bioinformatics/btv300 [PubMed]
55. Long J, Wang A, Bai Y, Lin J, Yang X, Wang D, Yang X, Jiang Y, Zhao H. Development and validation of a TP53-associated immune prognostic model for hepatocellular carcinoma. EBioMedicine. 2019; 42:363–74. https://doi.org/10.1016/j.ebiom.2019.03.022 [PubMed]
56. Gui J, Li H. Penalized cox regression analysis in the high-dimensional and low-sample size settings, with applications to microarray gene expression data. Bioinformatics. 2005; 21:3001–08. https://doi.org/10.1093/bioinformatics/bti422 [PubMed]
57. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010; 33:1–22. [PubMed]
58. Pencina MJ, D’Agostino RB Sr, D’Agostino RB Jr, Vasan RS. Evaluating the added predictive ability of a new marker: from area under the ROC curve to reclassification and beyond. Stat Med. 2008; 27:157–72. https://doi.org/10.1002/sim.2929 [PubMed]
59. Kiran M, Chatrath A, Tang X, Keenan DM, Dutta A. A prognostic signature for lower grade gliomas based on expression of long non-coding RNAs. Mol Neurobiol. 2019; 56:4786–98. https://doi.org/10.1007/s12035-018-1416-y [PubMed]

Research Paper Volume 13, Issue 2 pp 1929—1946