We retrieved the transcriptome and clinical profiles of training datasets from the The Cancer Genome Atlas (TCGA) database on the GDC data portal. The HTSeq-FPKM workflow was applied to compare all available bladder transitional epithelial carcinoma samples with normal bladder epithelium tissues. Then, we downloaded the validation datasets containing RNA sequencing data with expression profiles and clinical information from the Gene Expression Omnibus (GEO) database to verify the constructed model. We totally gathered 414 tumor samples and 19 normal tissues from the TCGA and 308 tumor samples from GEO (GSE32894). Then, 200 IRGs that drive, suppress, or mark inflammatory response were acquired from the gene set HALLMARK_INFLAMMATORY_RESPONSE (M5932) in the MsigDB molecular signature database, which are shown in Supplementary Table S1. The Human Protein Atlas database was used to observe immunohistochemistry (IHC) staining of genes with prognostic values.

Differentially expressed genes (DEGs) between BLCA and normal bladder tissues were classified and screened out by R package “limma” with a |log2 fold change (FC)| > 1 and a false discovery rate (FDR) < 0.05 in TCGA cohort. The information of sample sources (TCGA–BLCA) was included as covariates in the analysis. All screened out genes were included in the univariate Cox analysis to clarify the relationship between multiple risk factors and gene-related prognosis.

The LASSO algorithm was used to minimize the overfitting risk with “glmnet” R package, so as to adjust some regression coefficients return to zero. In LASSO regression, the independent variable was normalized expression matrix of candidate prognostic genes. The dependent variable was the OS of patients in the TCGA cohort. The 10-fold cross-validation was used to determine the penalty parameter of the prognostic model and was followed the minimum criteria. The risk scores of patients were weighted calculated according to the expression quantity of each IRG and its corresponding regression coefficient. The formula was as follows: R i s k s c o r e = ∑ i = 1 n C o e f ( x i ) ×   E x p ( x i )

In this formula, Coef (Xi) stood for the coefficient of each IRGs Xi, and Exp (Xi) stood for the expression levels of these genes. The risk scores were calculated and assigned to every individual. All involved individuals were assigned to the high- or low-risk groups, using the median risk score as a cutoff point. In line with the expressive abundance of genes in the signature, we used R package “Rtsne” and “ggplot2” to perform the PCA analysis and t-SNE analysis and thus to distinguish the distribution of individuals belong to different risk groups. The R package “survminer” was implemented to analyze the survival time of different cohort. The R package “survival” and “timeROC” were used to make time-independent receiver operating characteristic (ROC) curve to examine the efficacy of the prognostic model. Furthermore, Cox analyses were applied for selection of independent prognosis indicators.

Gene Ontology (GO, c5. go.v7.4. symbols) and Kyoto Encyclopedia of Genes and Genomes (KEGG, c2. cp.kegg.v7.4. symbols) analyses were performed with the “clusterProfiler” R package on the basis of signature genes using gene set enrichment analysis (GSEA) (version 4.1.0) software. Multiple GSEA analysis was completed to integrate the result of GO and KEGG analyses with “ggplot2” and “grid” R package.

The “GSVA” R package was used to perform a single-sample GSEA (ssGSEA) to clarify the immune scores of included samples. The correlation between immune scores and risk classification was calculated by operating R package “limma”. To visualize the relationship between filtration of immune cells and relevant immune signal pathways, box plots were drew with R package “ggpubr” and “reshape2”. The association between risk scores and immune infiltration subtype was tested by two-way ANOVA analysis and performed with “ggpubr” R package. The infiltration of immune cells was correlated with risk scores in different methods, including TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC.

The infiltration of immune cells and stromal cells in tissues of different samples were calculated and named as TME scores, including immune score and stromal score. The correlation between risk scores and TME scores was explored with spearman correlation analysis. The potential connection between risk groups and expression of common m6a epigenetic regulatory factors was calculated by R package “ggplot2” and “ggpubr”.

To examine whether the genes in signature were closely relevant to the drug sensitivity of chemotherapy, we utilized the NCI-60 database containing 60 different cancer cell lines from nine different types of tumors (Supplementary Table S2) and accessed through the CellMiner interface. Pearson correlation analysis was performed to investigate the correlation, which was visualized by C-map. Correlation analysis curves were built up with the efficacy of 218 drugs approved by Food and Drug Administration or in clinical trials (Supplementary Table S3).

R software (version 4.0.2) was applied to perform all statistical analyses and output visualization. Wilcoxon test was used to compare the DEGs between the pathological and normal bladder tissues. Univariate and multivariate Cox proportional hazard regression analyses were conducted to assess the association between risk score and prognosis. The ssGSEA immunity scores between the high- and low-risk groups were compared by Mann–Whitney test, while the p-values adjusted by the Benjamini–Hochberg method.

The study gathered a cohort of BLCA containing 414 tumor tissues and 19 normal tissues from TCGA database. Then, the validation group was constructed on the basis of GSE32894, which consists of 308 tumor tissues. Samples without survival time were eliminated and the amounts of the cohort were 407 in TCGA dataset and 224 in GSE32894. Some of the important clinical features of these patients are summarized in Table 1. The TCGA–BLCA cohort was set to be the trained group, whereas the validation group consolidated the GSE32894. Because we have extracted the expression values of 200 IRGs in patients with BLCA, 18 upregulated genes and 32 downregulated genes were authenticated (FDR < 0.05, log2 FC > 1). The differentially expressed IRGs were visualized via volcano plot (Figure 2A).

First, we excluded seven samples without complete clinical information, which was necessary for the subsequent analysis in the TCGA training group. Then, we used univariate Cox regression to screen out 10 IRGs as prognosis-related gene (Figure 2B). The construction of the gene–gene correlation map (Figure 2C) and protein–protein interaction (PPI) network of these gene-coding proteins (Figure 2D) revealed the combination and interaction among these genes. To prevent the risk of overfitting and consequent deviation, we utilized the LASSO–Cox regression analysis to optimize the prognostic signature associated with the expression of the 10 genes above, in an attempt to avoid overfitting the model and the consequent deviation. Then, a better consolidated model consisted of seven core genes was established by the minimum value of lambda (λ) (Figures 2E–G). The gene set was divided into protective genes (PTGER4 and RIPK2) and risky genes (LDLR, ROS1, MMP14, TNFAIP6, and MYC) on the basis of hazard ratio (HR). The details of these genes are summarized in Supplementary Table S4. Supplementary Figure S1 shows the relationship between the survival and screened IRGs. We further established a risk scoring model to facilitate the assessment and prediction of prognosis and long-term progression. The risk classification was demarcated by the median value of the whole population.

To verify the reliability and feasibility of the established prognostic signature, we drew a scatter plot to illustrate the distribution of risk scores (Figure 3A) and the correlation between risk scores and OS (Figure 3B) in the training group. The individuals in the low-risk group appeared to have apparently lower probability to encounter death ending than that in the high-risk group. We downloaded the GSE32894, which contained 308 tumor tissues, to reconfirm the conclusion via the same demarcation (84 samples were abandoned for lack of survival and status data) (Figures 3C,D). The clinical stratification analysis was used to explore whether the IRG signature could be widely applied to predict the survival and prognostic conditions in the age (≤60 and >60), gender (male and female), and TMN stage. The results indicated that the individuals in the low-risk group might perform a better prognosis than that in the high-risk group for every subgroup (p < 0.01, Supplementary Figures S2C-I). PCA and t-SNE plot illustrated that patients were distributed into two sections (Figures 3E–H). The survival curves showed a better survival in the low-risk cohort (p < 0.001) in the training group, whereas the validation group also support this conclusion (p < 0.01, Figures 3I,K). The ROC curves were used to evaluate the efficiency of risk scores in predicting the OS, and the area under the curve (AUC) of the seven-gene model was 0.666 in 1 year, 0.637 in 2 years, and 0.628 in 3 years. In the validation group, this value was 0.714 in 1 year, 0.691 in 2 years, and 0.681 in 3 years (Figures 3J, L). Consistently, the IHC staining also verified the expression level of signature genes in normal and tumor tissues. (Figure 4).

Furthermore, we analyzed the correlation between the risk scores and clinical characteristics of patients with BLCA by univariate and multivariate Cox regression analysis. The ability of the signature in predicting prognosis was independent. The univariate Cox regression identified the risk score [p < 0.001, HR = 3.227, 95% confidence interval (CI) = 2.081–5.004] and related clinical parameters including stage (p < 0.001, HR = 1.732, 95% CI = 1.427–2.101), grade (p = 0.141, HR = 2.855, 95% CI = 0.706–11.543), gender (p = 0.474, HR = 0.887, 95% CI = 0.639–1.232), and age (p < 0.001, HR = 1.033, 95% CI = 1.017–1.049) were, respectively, associated with OS in the TCGA training cohort (Figure 5A). Multivariate Cox regression analysis revealed that the risk score (p < 0.001, HR = 2.781, 95% CI = 1.755–4.408) was the independent prognostic factor of the OS (Figure 5B). Besides, similar to the results of the training group, the risk score (p = 0.022, HR = 4.821, 95% CI = 1.258–18.481) in the GSE32894 validation cohort was also proven to be the independent prognostic factors influencing the OS (Figures 5C,D; Table 2). The ROC curves were drawn to evaluate efficiency of the Cox regression analysis (Figures 5E,F). According to univariate and multivariate Cox regressions, we built a nomogram (Iasonos et al., 2008) to assess the probability of 1-, 2-, and 3-year OS of patients (Supplementary Figures S2A, 2B).

To further explore the biological functions and mechanistic pathways that were correlated with IRGs in BLCA, the GO and KEGG analysis were performed between the high-risk and low-risk groups in both TCGA and GEO cohorts. Interestingly, several inflammatory response–related molecular functions were enriched, such as extracellular matrix (ECM) organization, human immune response, and antimicrobial humoral response (Figures 6A,B). In addition, IRG set mainly participated in several metastasis-related biological processes including focal adhesion, ECM–receptor interaction, and PI3K-Akt signaling pathway (Figures 6C,D). The integration of the GO enrichment and KEGG pathway analysis came to a conclusion that the seven-gene signature was closely related to metastasis and invasion of tumor. The same analysis was performed on the basis of the GSE32894 dataset to verify the accuracy of the trained group. The result indicated that the IRG set was linked to epidermis development, cornification, and IL-17 signaling pathway (Supplementary Figure S3). Multiple GSEA analysis was utilized to explore the potential biological mechanisms of the inflammatory response–related signature involved in BLCA progression (Supplementary Figure S4).

To further investigate the immune feature in the high- and low-risk subgroups, the ssGSEA was performed to analyze the differences in the immune cell infiltrations and immune signal pathways. It was found that, in the high-risk group, the immune scores in DCs, macrophages, pDCs, Th1 cells, APC co-stimulation, and T cell co-inhibition were significantly higher. The score of Th2 cells and type Ⅱ interferon (IFN) response was significantly decreased in the high-risk group (Figures 7A–D). Subsequently, we explored the correlation of risk scores with immune infiltration, to confirm how risk score was associated with immune subtypes. Previous studies revealed that there are six primary types of immune infiltrations that were identified in human solid tumors. These subtypes, ranging from tumor progression to tumor suppression, respectively, were named as C1 (wound healing), C2 (INF-g dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-b dominant) (Tamborero et al., 2018). In the correlation plot of our study, almost all samples were concentrated from C1 to C4. Because no sample from TCGA–BLCA was matched to C5 and C6 immune subtypes, these two immune subtypes were excluded from the study. In addition, most individuals in the high-risk group aggregated in the subtype C1 and C2 (totally 93%) (Figures 7E,F), which indicated that C1 and C2 subtypes were correlated with poor prognosis. This result partially provided a brand-new research tendency and therapeutic guideline for immune infiltrated tumors. The infiltration of common immune cells was showed in heatmap on the basis of TCGA dataset (Supplementary Figure S5). As stromal cells play an essential role in the construction of the TME, particularly in MIBC, we conducted a correlation investigation of prognostic genes in the immune microenvironment. The IRG signature was found to be closely related to the stomal score of BLCA, indicating that the IRG sets were significantly positively characterized in relation to the TME formation (Figures 8B,C).

Then, we investigated the relevance of IRG signature to the content of immune checkpoint gene expression levels. In particular, PD-L1, as a key regulator in cancer immune evasion and surveillance, was an important indicator to guide individualized immunotherapy. As shown in Figures 8E,F, the risk score was significantly positively correlated with PD-L1 expression (r = 0.26, p < 0.001). The expression level of similar immune checkpoint genes was partly correlated with effect of immunotherapy. The result of immune checkpoint analysis revealed as expected that the expression of some famous targets was significantly higher in the high-risk group compared with the lower one (Figure 8D). Furthermore, we also investigated the relationship between m6a-related genes and risk scores. The m6A modifications were involved in almost all aspects of RNA metabolism, and transcriptome analysis of m6A-containing genes also revealed that the most relevant functional pathways were related to RNA metabolism (Sun et al., 2019). Previous studies revealed that m6a modification also greatly affected tumor progression. The correlation plots showed that the expression of m6a-related genes was elevated in the high-risk group (Figure 8A), illustrating that the signature was partly connected with epigenetic modification and took role in tumor progression and evolution.

To explore the relationship between the seven-gene signature and the BLCA drug resistance, we built the Scatter plots reflecting the expression distribution of the seven genes in some drug resistance datasets. Then, we investigated the expression of signature genes in NCI-60 cell line. The results revealed that high expression of MYC, MMP14, TNFAIP6, and PTGER4 were correlated with drug sensitivity of cancer cells to a multitude of chemotherapeutic drugs, especially MYC, which was highly correlated with Palbociclib, Carmustine, Ifosfamide, Parthenolide, etc. Conversely, increased expression of MMP14, RIPK2, LDLR, and PTGER4 was associated with increased drug resistance of cancer cells to Oxaliplatin, Lomustine, Decitabine, Dexrazoxane, trametinib, etc. (Figure 9).