Comprehensive Exploration of Tumor Microenvironment Modulation Based on the ESTIMATE Algorithm in Bladder Urothelial Carcinoma Microenvironment

Recently, the tumor microenvironment (TME) has been reported to be closely related to the tumor initiation, progression, and prognosis. Bladder urothelial carcinoma (BLCA), one of the most common subtypes of bladder cancer worldwide, has been associated with increased morbidity and mortality in the past decade. However, whether the TME status of BLCA contributes to the prediction of BLCA prognosis still remains uncertain. In this study, the ESTIMATE algorithms were used to estimate the division of immune and stromal components in 406 BLCA samples downloaded from The Cancer Genome Atlas database (TCGA). Based on the comparison between ESTIMATE scores, the differentially expressed genes (DEGs) were selected. Using the univariate Cox regression analysis, prognosis-related DEGs were further identified (p < 0.05). The LASSO regression analysis was then used to screen 11 genes that were highly related to the TME of BLCA to generate a novel prognostic gene signature. The following survival analyses showed that this signature could effectively predict the prognosis of BLCA. The clinical value of this signature was further verified in an external cohort obtained from the First Affiliated Hospital of Wenzhou Medical University (n = 120). Based on the stage-correlation analysis and differential expression analysis, IGF1 and MMP9 were identified as the hub genes in the signature. Additionally, using CIBERSORT algorithms, we found that both IGF1 and MMP9 were significantly associated with immune infiltration. Collectively, a novel TME-related prognostic signature contributes to accurately predict the prognosis of BLCA.

Introduction

Bladder urothelial carcinoma (BLCA), one of the most common subtypes of bladder cancer (BC) worldwide, has been associated with increased morbidity and mortality in the past decade (1, 2). Also, BLCA is the seventh most common cause of cancer-related mortality in the world. BLCA could be categorized into two distinctive types: muscle-invasive bladder cancer (MIBC) and nonmuscle-invasive bladder cancer (NMIBC), of which 25% represents MIBC and 75% represents NMIBC (3). Patients with NMIBC have a 50%–70% rate of relapse, with a 5-year survival rate >90%. However, after progressing to MIBC, the 5-year survival rate will drop to <50% (4). With the development of treatments and surgical techniques, patients with MIBC still have a poor prognosis with a 5-year survival rate of 5% (5, 6). Therefore, there is an urgent need to explore novel, effective biomarkers for the treatment of BLCA.

Recent studies have shown that the status of the tumor microenvironment (TME) is involved in the development of cancers (6). TME mainly contains tumor stromal cells and extracellular matrix (7). Increasing evidence has demonstrated that TME participates in the regulation of the invasiveness of tumor cells via chemotherapy, immune scape, etc. (8, 9). The types of nontumor components in TME mainly consist of immune cell and stromal cells, which play a key role in the diagnosis and prognosis of various human cancers (10). Immune cells have been reported to be associated with TME chemotaxis (11, 12). Moreover, stromal cells have been shown to be involved in tumor angiogenesis as well as extracellular matrix remodeling (13, 14). Notably, immune infiltration plays a crucial role in the modulation of TME status (15, 16).

Currently, to explore the novel TME-related prognostic signature for prognostic treatment of BLCA, there is an increase in interest in the evaluation of the stromal and immune components in the TME and tumor immune infiltration, especially in the comprehensive exploration of TME modulation. In this study, the stromal and immune scores were calculated via the ESTIMATE algorithm. The correlation analyses between ESTIMATE scores and clinical features revealed that TME status may be associated with the prognosis of BLCA. Based on the least absolute shrinkage and selection operator (LASSO) regression analysis and multivariate Cox regression analysis, a novel TME-related signature was generated to improve the prognosis prediction of BLCA. In addition, IGF1 and MMP9 were identified as the hub genes in the signature and used for the further immune infiltration analyses.

Materials and Methods

Data Collections

The cases of BLCA transcriptome data (FPKM normalized), including adjacent nontumorous samples and BLCA samples, were downloaded from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov). The relevant clinicopathological characteristics of BLCA patients were also downloaded from TCGA database (Table 1). According to each patient’s ID number, we matched the transcriptomic data and clinical features, and the data of the mismatched patients were removed. Finally, the gene expression profiles of 406 BLCA patients were obtained. An external cohort (named FAHWMU, n = 120) was obtained from the First Affiliated Hospital of Wenzhou Medical University (Wenzhou, China) to validate the clinical value of the prognostic signature. The collection of this cohort was reviewed and approved by the human research ethics committee of the First Affiliated Hospital of Wenzhou Medical University. The patients/participants provided their written informed consent to participate in this study.

TABLE 1

Table 1 Clinical characteristics of patients in TCGA database.

ESTIMATE Algorithms and Differential Expression Analyses

On the basis of the ESTIMATE algorithm, the relevant scores of infiltrating stromal cells and immune cells in BLCA patients were estimated (17). According to the median values of the immune and stromal scores, we divided the BLCA cases into two groups, the high group and the low group, respectively. We then obtained the heatmaps of the immune and the stromal scores using the R package named “heatmap.” The differential expression analyses were performed through the Wilcoxon test. The differentially expressed genes (DEGs) with p < 0.05 and | log fold change | >1.5 were filtered in the immune and stromal score groups. These genes with log fold change >1.5 were defined as upregulated genes, while genes with log fold change <−1.5 were defined as downregulated genes.

Enrichment Analyses

The DEGs were used for enrichment analyses from Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) databases using “ggplot2,” “org.Hs.eg.db,” “enrichplot,” and “clusterProfiler” packages. The DEGs with both p- and Q-values <0.05 were considered the significantly enriched genes.

Univariate Cox Regression Analysis and the Protein–Protein Interaction Network Construction

DEGs with p < 0.05 were screened as prognosis-related genes through the univariate Cox regression analysis. Using the Retrieval of Interacting Gene (STRING) database (https://www.string-db.org/), the protein–protein interaction (PPI) network was then obtained. The interactive relationships that were greater than 0.15 were used as the nodes. To better identify the hub genes, the Cytoscape software (version 3.8.2) (https://cytoscape.org/) was used to reconstruct the PPI network. The Gene Set Enrichment Analysis (GSEA) software (version 4.1.0) (http://www.gsea-msigdb.org/gsea/) was used to analyze all tumor cases through C7.all.v7.4.symbols.gmt gene sets.

Generation of the Prognostic Signature

The LASSO regression analysis was performed among 45 DEGs regarding overall survival (OS) to remove the ones that were overfitted to the signature. The multivariate Cox regression analysis was then used to identify 11 hub genes related to the BLCA microenvironment, and their coefficients were obtained to calculate the risk score. The risk score was calculated according to the following computational formula: risk score = Ʃ (β_i * Exp_i) (i = 11). In addition, according to the median values of risk scores, we categorized the BLCA patients into the high-risk and the low-risk groups. The “survival” and “survival ROC” packages were used to generate the receiver operating characteristic curve (ROC) to predict the prognostic value of the signature. The log rank test was used to perform the Kaplan–Meier survival analysis between the high-risk and low-risk groups.

Identification of the Hub Genes

The “ggpubr” and “limma” R packages were used to analyze the correlations between TNM stage and the genes of the signature. Then, the hub DEGs between adjacent nontumorous samples and BLCA samples were screened through differential expression analysis. Lastly, to determine the hub genes of the prognostic signature, the “Venndiagram” R package was used to identify the common genes between stage-correlation analyses and differential analyses.

Immune Infiltration Analyses

In BLCA patients, the contents of immune-infiltrating cells were obtained via the CIBERSORT algorithm. Only samples with a CIBERSORT p < 0.05 were included for the next analysis. Subsequently, “ggExtra,” “vioplot,” “ggpubr,” and “gglot2” R packages were applied to perform the stage-correlation analyses. In the analysis of Pearson coefficient tests and Wilcoxon rank-sum, p < 0.05 was set as the criterion for a significant difference. In addition, immune correlation analysis was performed using the “limma” R package and Person’s correlation coefficients.

Results

The Correlation Between Clinicopathological Characteristics and the TME of BLCA Was Significant

The overall workflow of this study is shown in Figure 1. Recently, several studies have shown that a higher immune or stromal score indicates more immune or stromal components in the TME (18). Using the ESTIMATE algorithm, it was found that the ESTIMATE, the stromal, and the immune scores were significantly associated with gender and grade, especially for the T2 and T3 histologically graded disease (Figures 2A–I). The p-values of gender were 0.036, 0.034, and 0.028, respectively, while the p-values of grade were 0.00028, 7.4e−07, and 5.3e−06, respectively. In addition, the p-values between T2 and T3 were 0.067, 3.4e−07, and 0.00022, respectively. The stromal scores were also positively associated with N0 and N2 of histological grades (p = 0.0018). However, no significant relationships between the scores and other clinicopathological features were found (Supplementary Figure S1). These correlation analyses suggest that the stromal and immune scores may be associated with the progression of BLCA.

FIGURE 1

Figure 1 The overall workflow of this study.

FIGURE 2

Figure 2 Significant correlation between clinicopathological characteristics and scores. (A–C) Correlation between the scores and grade. (A) Distribution of ImmuneScore in grade (p = 0.00028). (B) Distribution of StromalScore in grade (p = 7.4e−07). (C) Distribution of ESTIMATEScore in grade (p = 5.3e−06). (D–F) Correlation between scores and T classification. (D) Distribution of ImmuneScore in T classification. (E) Distribution of StromalScore in T classification. (F) Distribution of ESTIMATEScore in T classification. (G–I) Correlation between scores and gender. (G) Distribution of ImmuneScore in gender (p = 0.036). (H) Distribution of StromalScore in gender (p = 0.034). (I) Distribution of ESTIMATEScore in gender (p = 0.028).

A Total of 1,097 Genes Were Identified as DEGs

Based on the median values of the immune score as well as stromal score, all BLCA patients were divided into two groups. As shown in the heatmaps, DEGs were identified by the differential analyses (Figures 3A, B). On the basis of the median values of immune score, there were 1,662 DEGs including 384 downregulated genes and 1,278 upregulated genes. Similarly, based on the stromal score, 1,613 DEGs including downregulated 199 and upregulated 1,414 were identified. The common DEGs between the immune score group and the stromal score group were shown in the Venn diagram (Figures 3C, D). Finally, a total of 1,097 common DEGs were identified and used for the next studies.

FIGURE 3

Figure 3 Identification of 1,027 DEGs between ImmuneScore and StromalScore. (A) Heatmap for DEGs in ImmuneScore, only showing the top 60 significantly DEGs (p < 0.05, fold change > 1.5). (B) Heatmap for DEGs in StromalScore, only showing the top 60 significant DEGs (p < 0.05, fold change > 1.5). (C) Venn plots showing 109 DEGs that were commonly upregulated both in ImmuneScore and StromalScore. (D) Venn plots showing 988 DEGs that were commonly downregulated both in ImmuneScore and StromalScore.

Immune-Related Pathways Were Enriched in the GO and KEGG Enrichment Analyses

The GO and KEGG pathway enrichment analyses were performed to select the potential molecular pathways among 1,097 DEGs. GO analysis showed that the DEGs were mainly enriched in humoral immune response mediated by circulating immunoglobulins, lymphocyte-mediated immunity, and immune response-activating cell surface receptor signaling pathway. These pathways were strongly associated with immune system (Figure 4A). Moreover, results of KEGG analysis indicated that the DEGs were significantly associated with cytokine–cytokine receptor interactions and immune pathways (B-cell receptor signaling pathway, T-cell receptor signaling pathway, and intestinal immune network for IgA production) (Figure 4B). Our results suggest that immune-related pathways may participate in the progression of BLCA.

FIGURE 4

Figure 4 GO and KEGG enrichment analyses. (A) GO enrichment analysis for 1,097 DEGs, only both p and q < 0.05 were considered significant. (B) KEGG pathways enrichment analysis for 1,097 DEGs, only both p and q < 0.05 were considered significant.

Identification of the Hub Genes

Next, the 1,097 DEGs were exploited into univariate Cox regression analysis. Thus, 45 prognostic-related DEGs were identified (p < 0.05) (Figure 5). Using the SRING tool and Cytoscape software, the PPI network was constructed to explore their underlying interplays (Figure 6A). During PPI network, many genes with too few nodes (nodes < 2) were deleted. As shown in Figure 6B, DCN, IGF1, and MMP9 were identified as relatively remarkable nodes. Furthermore, GSEA enrichment analysis revealed that 3 immune-related pathways including the CD4 T-cell pathway, MEM CD4 T-cell pathway, and native B-cell pathway, may influence the progression of BLCA (Figures 6C–E).

FIGURE 5

Figure 5 The univariate COX regression analysis of the DEGs. Only 45 DEGs were significantly correlated with prognosis (p < 0.05).

FIGURE 6

Figure 6 GSEA enrichment analyses of the hub genes. (A) The PPI network of 45 PRDEGs (coefficients > 0.15). (B) DCN, MMP9, and IGF1 had the most nodes in the PPI network. (C–E) The GSEA enrichment analyses of IGF1 (C), MMP9 (D), and IGF1 (E).

Construction and Verification of the Prognostic Signature

The LASSO analysis was performed to avoid overfit of the signature (Figures 7A, B). Subsequently, 11 TME-related genes were identified by the multivariate Cox analysis (Figure 7C). The risk formula was calculated as:

FIGURE 7

Figure 7 Generation of the prognostic risk signature. (A, B) The LASSO regression analysis of the 45 PRDEGs. (C) The Forest plot showing eight hub genes selected by the multivariate regression analysis.

Risk score = −0.414 * ADGRG5 + 0.001 * MMP9 + 0.222 * IGF1 − 0.526 * AC015911.3 + 0.726 * AL133415.1 − 0.753 * LINC02362 − 0.123 * AC078880.3 − 0.050 * AC012307.1 − 0.050 * FAHD2P1 − 0.034 * LINC00930 − 1.220 * SPINK14.

A total of 406 BLCA patients were then divided into the low-risk group (n = 203) and high-risk group (n = 203) according to the median risk score (Figure 8C). The distribution of survival status of the low-risk group as well as the high-risk group is shown in Figure 8D. The area under the curve (AUC) value of ROC for 3-year OS prediction reached 0.749, indicating that the signature had a good prediction for BLCA prognosis (Figure 8B). In addition, it was found that the OS of BLCA patients in the high-risk group was obviously worse than that in the low-risk group (Figure 8A). The prediction accuracy of the signature in BLCA prognosis was validated in the FAHWMU cohort (n = 120). Similarly, BLCA patients at low risk had better OS compared with those at high risk (Supplementary Figure S2A). The AUC was 0.793 in the 1st year, 0.752 in the 2nd year, and 0.810 in the 3rd year, respectively (Supplementary Figure S2B).

FIGURE 8

Figure 8 Survival analyses of the signature. (A) The Kaplan-Meier survival curve analysis showing the significant difference between the high-risk and the low-risk groups (p < 0.05). (B) Time-dependent ROC curve analysis showing the significant prognostic value (all AUC value > 0.70). (C) Risk score distribution of the patients. (D) Survival status scatter plots showing poorer prognosis of the high-risk group compared with the low-risk group.

IGF1 and MMP9 Were Identified as the Hub Genes in the Prognostic Signature

As indicated by the clinical clustering heatmap, this TME-related prognostic signature was significantly associated with TNM stage and grade (Figure 9A, p < 0.001). The correlation analysis showed that IGF1 and AL133415.1 were positively correlated with TNM stage. By contrast, a negative correlation was found between TNM stage and SPINK14 as well as MMP9, AC078880.3, LINC00930, and AC012307.1 (Figures 9B–H, p < 0.01). No significant correlations were found between the TNM stage and other genes, including AC015911.3, ADGRG5, FAHD2P1, and LINC02362 (Supplementary Figure S3, p > 0.01). Moreover, in comparison with the adjacent nontumorous samples, MMP9 and FAHD2P1 were enhanced in the BLCA samples, whereas IGF1 was reduced (Figures 10A–C, p < 0.01). No significant difference was found in other genes between the BLCA and nontumorous samples (Supplementary Figure S4, p > 0.01). As shown in the Venn diagram, two common genes (IGF1 and MMP9) were selected through the combination of stage-correlation analyses and differential analyses (p < 0.01) (Figure 10D). Thus, IGF1 and MMP9 were identified as the hub genes and selected for the next studies.

FIGURE 9

Figure 9 Stage correlation of the risk score and 11 hub genes in the signature. (A) Multifactorial heatmap showing certain clinical features with a significant correlation with the risk score. (^*p < 0.05, ^**p < 0.01, ^***p < 0.001) (B–H). The expression of IGF1 and AL133415.1 were positively correlated with TNM stages, while the expression of SPINK14, MMP9, AC078880.3, LINC00930, and AC012307.1 were negatively associated with the TNM stages.

FIGURE 10

Figure 10 IGF1 and MMP9 were identified as the hub genes in the signature. (A, B) The expression of MMP9 and FAHD2P1 was significantly higher in the tumor samples than in the normal samples (p < 0.01). (C) The expression of IGF1 was significantly lower in the tumor samples compared with the normal samples (p < 0.001). (D) The Venn plot showing IGF1 and MMP9 as the common genes between DEGs and stage-correlated genes.

IGF1 and MMP9 Were Involved in the Immune Activity of TME in BLCA

The fraction of 22 tumor-infiltrated immune cells in the TCGA cohort was calculated through CIBERSOFT algorithm (Figures 11A, B, p < 0.05). The intersection of differential analyses revealed that M2 macrophages were positively correlated with the expressions of IGF1 and MMP9. Other immune cells such as B-cell native, were negatively correlated with the expression of MMP9 (p = 0.007) and positively associated with the expression of IGF1 (p = 0.003) (Figures 12A, 13A). Furthermore, the correlation analyses revealed that several types of immune cells were obviously associated with the expressions of IGF1 and MMP9. With the increase in the expression of IGF1, native B-cell count was increased (R = 0.35, p < 0.05) (Figure 12B). Moreover, MMP9 expression was negatively correlated with the levels of native B cell (R = −0.18, p < 0.05) (Figure 13B). Finally, as shown in the Venn diagram, the common immune cells were selected by the differential analyses and immune cell correlation analyses (Figures 12C, 13C). All these results suggest that IGF1 and MMP9 are closely related to immune infiltrates and may be involved in the immune activity of TME in BLCA.

FIGURE 11

Figure 11 The tumor-infiltrated immune cell profile in BLCA samples and correlation analysis. (A) Barplot showing the proportion of 21 types of tumor-infiltrating immune cells in BLCA samples. (B) Heatmap showing the correlation between the 22 types of tumor-infiltrating immune cells.

FIGURE 12

Figure 12 Significant correlation between the expression of IGF1 and the fractions of tumor-infiltrating immune cells. (A) Violin plot showing the comparison between 22 types of tumor-infiltrating immune cells with IGF1 expression. (B) The correlation tests showing the correlation between 10 types of tumor-infiltrating immune cell proportion with the IGF1 expression (all p < 0.05). (C) Venn plot showing the collective TICs between the difference and correlation tests.

FIGURE 13

Figure 13 Significant correlation between the expression of MMP9 and the fractions of tumor-infiltrating immune cells. (A) Violin plot showing the comparisons between 22 types of tumor-infiltrating immune cells with MMP9 expression. (B) The correlation tests showing the correlation between 10 types of tumor-infiltrating immune cells proportion with the MMP9 expression (all p < 0.05). (C) Venn plot showing the collective TICs between the difference and correlation tests.

Discussion

Recently, the remodeling of the TME has been shown to play a crucial role in the initiation and progression of various human cancers (19). Increasing studies have demonstrated that TME modulation is associated with BLCA progression (20). However, whether the TME status of BLCA contributes to the prediction of BLCA prognosis remains unknown. Recent studies have been reported that the immune and stromal scores in the TME could be calculated by ESTIMATE algorithms, which contributes to provide a novel insight into the development of immune treatment (21). Currently, immunotherapy for TME has been shown to improve BLCA treatment (2). Therefore, more novel candidates are needed for BLCA immunotherapy (22).

In this study, we used the data from the TCGA database to calculate the ESTIMATE, stromal, and immune scores. Our results showed that immune scores and stromal scores were significantly higher in BLCA, associated with lower TNM stages as well as higher stage T levels. The GO and KEGG enrichment analyses indicated that the DEGs were mainly enriched in immune-related pathways. The GSEA analysis confirmed that three hub genes including DCN, IGF1, and MMP9, had a strong correlation with the immune pathway. Moreover, survival analyses preliminarily revealed that the signature could accurately predict the prognosis of BLCA patients. With the aid of the FAHWMU cohort, the clinical value of our prognostic signature was further confirmed.

Recently, TME-related signatures have been identified and validated in multiple human cancers. Huang et al. identified a novel long noncoding RNA signature to improve the prognosis and immunotherapy response of patients with hepatocellular carcinoma (23). Using multiomics analysis, Shi et al. determined a TME-related signature to identify distinct prognostic patterns in osteosarcoma (24). Similarly, the signature correlated with TME was also constructed in the renal cancer to improve immunotherapeutic and antitumor approaches (25). Compared with previous studies, the novel 11 TME-related gene prognostic signature determined in this study has many advantages. The prognostic signature showed a good prognostic value, which was further confirmed in an external BLCA cohort. In addition, the risk scores of the signature were found to be an independent prognostic factor for BLCA treatment. Moreover, as the hub genes in the signature, IGF1 and MMP9 were significantly correlated with immune infiltration. Taken together, our results suggest that this novel TME-related prognostic signature may provide a new theoretical basis for the prognostic assessment of patients with BLCA.

Herein, IGF1 and MMP9 were confirmed as the independent prognostic indicators for BLCA. IGF1, a polypeptide that contributes to tumor progression, has been reported to be involved in TME status (26). Young et al. found that the low expression of IGF1 could induce the aging of hematopoietic stem cells (26). Shi et al. found that IGF1 is associated with the expressions of immune cells, such as NK cells (27). In this study, we found that IGF1 may contribute to an increase in the levels of native B cell and M2 macrophages. MMP9, as a classical immune-related gene, has been demonstrated to play a vital role in the immune treatment of BLCA (28). MMP9 contributes to the promotion of the neutrophil-mediated inflammation, indicating its roles in immune infiltrates (29). Consistent with the previous studies, we found that MMP9 plays a crucial role in the regulation in the activities of immune cells such as native B cells and M0 macrophages, etc. Furthermore, using the STRING database, the potential functional link between TME-related genes and IGF1 as well as MMP9 was obtained (Supplementary Figure S5). Combined with these, IGF1 and MMP9, the hub genes of our prognostic signature, may participate in the progression of BLCA via the immune activity of TME.

However, there are many limitations in this study. The prognostic value of this signature still needs to be verified in clinical samples with large numbers. In addition, the underlying mechanisms of the 11 TME-related hub genes in immune infiltrates of BLCA should be explored in the future.

In conclusion, we disclose a novel 11 TME-related genes prognostic signature for BLCA, which may provide a new theoretical basis for the prognosis of BLCA patients.

Data Availability Statement

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

Ethics Statement

The studies involving human participants were reviewed and approved by the Human Research Ethics Committee of The First Affiliated Hospital of Wenzhou Medical University. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author Contributions

JC and CL designed the study and analyzed the data. JC, BL, YZ, KZ, and RZ revised the images. JC, BC, and YJ performed the literature search and collected data for the manuscript. JZ and YL revised the manuscript. All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

The project was supported by Wenzhou Medical University Basic Scientific Research (No. KYYW201904).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

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

Supplementary Material

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

Supplementary Figure 1 | The correlations of other clinicopathological characteristics with scores. (A) The distribution of score in age. (B) The distribution of scores in M classification. (C) The distribution of scores in M classification. (D) The distribution of score in TNM stage.

Supplementary Figure 2 | Survival analysis of the external cohort. (A) Kaplan-Meier survival curve of high- and low-risk patients. (B) ROC curve of the signature at 1^st, 2^nd and 3^rd year.

Supplementary Figure 3 | The correlations between the other genes and TNM stage were not significant (P > 0.01).

Supplementary Figure 4 | The expression levels of the other genes in the signature between BLCA samples and normal samples (P > 0.01).

Supplementary Figure 5 | The network regarding the potential functional link between TME-related genes and IGF1 as well as MMP9.

Abbreviations

BLCA, bladder urothelial carcinoma; TME, tumor microenvironment; TCGA, The Cancer Genome Atlas; LASSO, least absolute shrinkage and selection operator; DEGs, differentially expressed genes; OS, overall survival; AUC, area under the curve; ROC, receiver operating characteristic; PPI, protein–protein interaction; GSEA, Gene Set Enrichment Analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Wang H, Wu X, Chen Y. Stromal-Immune Score-Based Gene Signature: A Prognosis Stratification Tool in Gastric Cancer. Front Oncol (2019) 9:1212. doi: 10.3389/fonc.2019.01212

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Babjuk M, Böhle A, Burger M, Capoun O, Cohen D, Compérat E, et al. EAU Guidelines on Non-Muscle-Invasive Urothelial Carcinoma of the Bladder: Update 2016. Eur Urol (2017) 71(3):447–61. doi: 10.1016/j.eururo.2016.05.041

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Rubben H, Lutzeyer W, Fischer N, Deutz F, Lagrange W, Giani G. Natural History and Treatment of Low and High Risk Superficial Bladder Tumors. J Urol (1988) 139(2):283–5. doi: 10.1016/s0022-5347(17)42387-1

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Babjuk M, Oosterlinck W, Sylvester R, Kaasinen E, Böhle A, Palou-Redorta J, et al. EAU Guidelines on non-Muscle-Invasive Urothelial Carcinoma of the Bladder, the 2011 Update. Actas Urologicas Espanolas (2012) 36(7):389–402. doi: 10.1016/j.acuro.2011.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Bejarano L, Jordao MJC, Joyce JA. Therapeutic Targeting of the Tumor Microenvironment. Cancer Discov (2021) 11(4):933–59. doi: 10.1158/2159-8290.CD-20-1808

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Thomas D, Radhakrishnan P. Pancreatic Stellate Cells: The Key Orchestrator of The Pancreatic Tumor Microenvironment. Adv Exp Med Biol (2020) 1234:57–70. doi: 10.1007/978-3-030-37184-5_5

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Cheng H, Lee J, Wahli W, Tan N. Exploiting Vulnerabilities of Cancer by Targeting Nuclear Receptors of Stromal Cells in Tumor Microenvironment. Mol Cancer (2019) 18(1):51. doi: 10.1186/s12943-019-0971-9

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Park J, van Wyk H, Roxburgh C, Horgan P, Edwards J, McMillan D. Tumour Invasiveness, the Local and Systemic Environment and the Basis of Staging Systems in Colorectal Cancer. Br J Cancer (2017) 116(11):1444–50. doi: 10.1038/bjc.2017.108

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Wörmann S, Diakopoulos K, Lesina M, Algül H. The Immune Network in Pancreatic Cancer Development and Progression. Oncogene (2014) 33(23):2956–67. doi: 10.1038/onc.2013.257

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Xiang X, Wang J, Lu D, Xu X. Targeting Tumor-Associated Macrophages to Synergize Tumor Immunotherapy. Signal Transduction Targeted Ther (2021) 6(1):75. doi: 10.1038/s41392-021-00484-9

CrossRef Full Text | Google Scholar

12. Bunse M, Pfeilschifter J, Bluhm J, Zschummel M, Joedicke J, Wirges A, et al. CXCR5 CAR-T Cells Simultaneously Target B Cell Non-Hodgkin's Lymphoma and Tumor-Supportive Follicular T Helper Cells. Nat Commun (2021) 12(1):240. doi: 10.1038/s41467-020-20488-3

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Tian C, Huang Y, Clauser K, Rickelt S, Lau A, Carr S, et al. Suppression of Pancreatic Ductal Adenocarcinoma Growth and Metastasis by Fibrillar Collagens Produced Selectively by Tumor Cells. Nat Commun (2021) 12(1):2328. doi: 10.1038/s41467-021-22490-9

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Pozzi S, Scomparin A, Israeli Dangoor S, Rodriguez Ajamil D, Ofek P, Neufeld L, et al. Meet Me Halfway: Are In Vitro 3D Cancer Models on the Way to Replace In Vivo Models for Nanomedicine Development? Advanced Drug Deliv Rev (2021) 175:113760. doi: 10.1016/j.addr.2021.04.001

CrossRef Full Text | Google Scholar

15. Qiao M, Jiang T, Liu X, Mao S, Zhou F, Li X, et al. Immune Checkpoint Inhibitors in EGFR-Mutated NSCLC: Dusk or Dawn? J Thorac Oncol Off Publ Int Assoc Study Lung Cancer (2021) 16(8):1267–88. doi: 10.1016/j.jtho.2021.04.003

CrossRef Full Text | Google Scholar

16. Jia L, Gao Y, Zhou T, Zhao X, Hu H, Chen D, et al. Enhanced Response to PD-L1 Silencing by Modulation of TME via Balancing Glucose Metabolism and Robust Co-Delivery of siRNA/Resveratrol With Dual-Responsive Polyplexes. Biomaterials (2021) 271:120711. doi: 10.1016/j.biomaterials.2021.120711

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Liu W, Ye H, Liu YF, Xu CQ, Zhong YX, Tian T, et al. Transcriptome-Derived Stromal and Immune Scores Infer Clinical Outcomes of Patients With Cancer. Oncol Lett (2018) 15(4):4351–7. doi: 10.3892/ol.2018.7855

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Liu D, Yang X, Wu X. Tumor Immune Microenvironment Characterization Identifies Prognosis and Immunotherapy-Related Gene Signatures in Melanoma. Front Immunol (2021) 12:663495. doi: 10.3389/fimmu.2021.663495

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zhang H, Song J, Dong J, Liu Z, Lin L, Wang B, et al. Tumor Microenvironment Analysis Identified Subtypes Associated With the Prognosis and the Tumor Response to Immunotherapy in Bladder Cancer. Front Genet (2021) 12:551605. doi: 10.3389/fgene.2021.551605

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Liu Y, Wu Y, Zhang P, Xu C, Liu Z, He C, et al. CXCL12 and CD3E as Indicators for Tumor Microenvironment Modulation in Bladder Cancer and Their Correlations With Immune Infiltration and Molecular Subtypes. Front Oncol (2021) 11:636870. doi: 10.3389/fonc.2021.636870

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhang LH, Li LQ, Zhan YH, Zhu ZW, Zhang XP. Identification of an IRGP Signature to Predict Prognosis and Immunotherapeutic Efficiency in Bladder Cancer. Front Mol Biosci (2021) 8:607090

Google Scholar

23. Huang S, Zhang J, Lai X, Zhuang L, Wu J. Identification of Novel Tumor Microenvironment-Related Long Noncoding RNAs to Determine the Prognosis and Response to Immunotherapy of Hepatocellular Carcinoma Patients. Front Mol Biosci (2021) 8:781307. doi: 10.3389/fmolb.2021.781307

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Shi D, Mu S, Pu F, Liu J, Zhong B, Hu B, et al. Integrative Analysis of Immune Related Multi-Omics Profiles Identifies Distinct Prognosis and Tumor Microenvironment Patterns in Osteosarcoma. Mol Oncol (2021) 10.1002. doi: 10.1002/1878-0261.13160

CrossRef Full Text | Google Scholar

25. Feng T, Zhao J, Wei D, Guo P, Yang X, Li Q, et al. Immunogenomic Analyses of the Prognostic Predictive Model for Patients With Renal Cancer. Front Immunol (2021) 12:762120. doi: 10.3389/fimmu.2021.762120

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Young K, Eudy E, Bell R, Loberg M, Stearns T, Sharma D, et al. Decline in IGF1 in the Bone Marrow Microenvironment Initiates Hematopoietic Stem Cell Aging. Cell Stem Cell (2021) 28(8):1473–82. doi: 10.1016/j.stem.2021.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Shi JW, Yang HL, Lai ZZ, Shen HH, Qin XY, Qiu XM, et al. WISP2/IGF1 Promotes the Survival of DSCs and Impairs the Cytotoxicity of Decidual NK Cells. Reproduction (2021) 161(4):425–36. doi: 10.1530/REP-20-0658

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Qu G, Liu Z, Yang G, Xu Y, Xiang M, Tang C. Development of a Prognostic Index and Screening of Prognosis Related Genes Based on an Immunogenomic Landscape Analysis of Bladder Cancer. Aging (Albany NY) (2021) 13(8):12099–112. doi: 10.18632/aging.202917

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang Q, Meng Y, Wang K, Zhang X, Chen W, Sheng J, et al. Inflammation and Antiviral Immune Response Associated With Severe Progression of COVID-19. Front Immunol (2021) 12:631226. doi: 10.3389/fimmu.2021.631226

PubMed Abstract | CrossRef Full Text | Google Scholar