{'fraction_non_alphanumeric': 0.0377906976744186, 'fraction_numerical': 0.023147141925668333, 'mean_word_length': 5.119534497090607, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification and Quantification of Iron Metabolism Landscape on Therapy and Prognosis in Bladder Cancer

The morbidity of bladder cancer (BLCA) is high and has gradually elevated in recent years. BLCA is also characterized by high recurrence and high invasiveness. Due to the drug resistance and lack of effective prognostic indicators, the prognosis of patients with BLCA is greatly affected. Iron metabolism is considered to be a pivot of tumor occurrence, progression, and tumor microenvironment (TME) in tumors, but there is little research in BLCA. Herein, we used univariate COX regression analysis to screen 95 prognosis-related iron metabolism-related genes (IMRGs) according to transcription RNA sequencing and prognosis information of the Cancer Genome Atlas (TCGA) database. TCGA-BLCA cohort was clustered into four distinct iron metabolism patterns (C1, C2, C3, and C4) by the non-negative matrix factorization (NMF) algorithm. Survival analysis showed that C1 and C3 patterns had a better prognosis. Gene set variant analysis (GSVA) revealed that C2 and C4 patterns were mostly enriched in carcinogenic and immune activation pathways. ESTIMATE and single sample gene set enrichment analysis (ssGSEA) also confirmed the level of immune cell infiltration in C2 and C4 patterns was significantly elevated. Moreover, the immune checkpoint genes in C2 and C4 patterns were observably overexpressed. Studies on somatic mutations showed that the tumor mutation burden (TMB) of C1 and C4 patterns was the lowest. Chemotherapy response assessment revealed that C2 pattern was the most sensitive to chemotherapy, while C3 pattern was the most insensitive. Then we established the IMRG prognosis signature (IMRGscore) by the least absolute shrinkage and selection operator (LASSO), including 13 IMRGs (TCIRG1, CTSE, ATP6V0A1, CYP2C8, RNF19A, CYP4Z1, YPEL5, PLOD1, BMP6, CAST, SCD, IFNG, and ASIC3). We confirmed IMRGscore could be utilized as an independent prognostic indicator. Therefore, validation and quantification of iron metabolism landscapes will help us comprehend the formation of the BLCA immunosuppressive microenvironment, guide the selection of chemotherapeutic drugs and immunotherapy, and predict the prognosis of patients.


INTRODUCTION
Bladder cancer (BLCA) is one of the most familiar malignant tumors in the urinary system, with about 81400 new cases and 17900 deaths in the United States in 2020 (Siegel et al., 2020). Approximately 75% of BLCA was found to be non-muscle invasive bladder cancer (NMIBC), which was characterized by a high recurrence rate (45% 5-year recurrence rate) (Berdik, 2017;Babjuk et al., 2019). Transurethral resection of bladder tumor (TURBT), chemotherapy, BCG vaccine, radiotherapy, and radical cystectomy are the main treatments for BLCA patients (Berdik, 2017). Chemotherapy and immunotherapy are also important strategies for conservative treatment of BLCA (Yin et al., 2016;Rouanne et al., 2018). However, some patients are not sensitive to these drug therapies. Due to the high recurrence rate, high metastatic risk, and patient's dissatisfaction with the treatment effect, it is of great significance to identify and quantify some molecular landscapes that have impacts on the choice of drug treatment, and explore a novel indicator that predicts the prognosis of BLCA patients.
Iron is a vital trace element for cell proliferation and growth in the human body (Torti and Torti, 2013). In cancer, the absorbability, effusion, storage, and regulation of iron are entirely disturbed, which indicates that the reprogramming of iron metabolism would induce the dysregulation in tumor cells division and survival (Andrews, 2008;Manz et al., 2016;Wang et al., 2018;Jung et al., 2019). Iron plays a dual role in cancer (Thévenod, 2018). Epidemiological investigations revealed that excess iron was a hazard factor of carcinogenesis (Stevens et al., 1994;Wu et al., 2004;Fonseca-Nunes et al., 2014). The accumulation of iron supports tumor worsening in proliferation, metabolism, and metastasis (Torti et al., 2018). Cancer cells exhibit a phenotype search for iron through disordering regulation of iron-binding proteins (Dufès et al., 2013;Bialasek et al., 2019). On the other hand, iron dependence of cancer cells affects many cell death modes, including ferroptosis, a form of irondependent cell death (Mou et al., 2019;Battaglia et al., 2020). Inducing ferroptosis of cancer cells has become a new hotspot in the research and development of cancer treatment (Hassannia et al., 2019;Xu et al., 2019).
There has been little research on iron metabolism in BLCA. The study was conducted to confirm whether iron metabolism had an effect on the molecular microenvironment of BLCA, as well as its ability to predict the clinical prognosis. We first clustered the TCGA-BLCA cohort into different iron metabolic patterns on the basis of the expression of iron metabolism-related genes (IMRGs). Then the survival prognosis, GSVA analysis, tumor immune microenvironment (TIME), somatic mutations chemotherapy, and immunotherapy response among different patterns were analyzed. Eventually, we established a prognostic signature associated with iron metabolism and confirmed that it is an effective independent predictor in BLCA patients.

NMF Clustering for Iron Metabolism Patterns
We matched IMRGs' RNAseq data and overall survival (OS) information of the TCGA-BLCA dataset. A univariate Cox  regression analysis was carried out to determine the prognosisrelated IMRGs in BLCA, with a screening criterion of p < 0.05. Non-negative matrix factorization (NMF) is applied to determine distinct iron metabolism-related patterns with the help of the "NMF" R package. NMF algorithm decomposes the original matrix into two non-negative matrices to identify the potential features in the gene expression profile (Brunet et al., 2004). Repeat the deposition and aggregate the results to obtain consistent clustering. According to the cophenetic coefficient, contour, and sample size, k = 4 was determined as the best cluster number.

GSVA
Gene set variant analysis (GSVA), is a nonparametric, unsupervised algorithm. GSVA transforms the isolate gene expression matrix to an expression matrix of particular gene sets as features. The algorithm is implemented based on the "GSVA" R package. The difference of expression matrix after transformation was analyzed by the "limma" package to find the difference of enriched functions among different iron metabolic patterns.

Evaluation of Tumor Immune Microenvironment
In order to assess the TIME status of BLCA, we used single sample gene set enrichment analysis (ssGSEA), ESTIMATE, and CIBERSORT in R. ssGSEA investigated congenital and adaptive immune cells as well as a variety of immune-related functions. The Normalized Enrichment Score (NES) was to embody the relative amount of each TIME infiltration unit in patients. ESTIMATE predicted the level of infiltrating matrix and immune cells by calculating stromal and immune scores, comprehensively obtained the ESTIMATE score for evaluating the TIME. We also evaluated the relative fraction of 22 tumorinfiltrating immune cells (TIICs), including B cells, T cells, natural killer (NK) cells, macrophages, dendritic cells (DCs), eosinophils, neutrophils, and so on in each cancer sample with CIBERSORT algorithm. p < .05 was the threshold of a credible sample for estimating the proportion of immune cells.

Evaluation of Immunotherapy and Chemotherapy on Iron Metabolism Patterns
Tumor Immune Dysfunction and Exclusion (TIDE) score is a computational framework developed based on the analysis and modeling of characteristic genes for T cell exclusion and T cell dysfunction in immunosuppression at high levels of cytotoxic T lymphocytes (CTL) (Jiang et al., 2018). We applied four indicators to predict the efficacy of immunotherapy, including exclusion score, dysfunction score, microsatellite instability (MSI), and TIDE. The chemotherapeutic response of BLCA patients was evaluated by Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerRxgene.org). Six chemotherapeutic drugs in BLCA treatment were chosen, including Gemcitabine, Cisplatin, Docetaxel, Mitomycin-C, Doxorubicin, and Paclitaxel. The ridge regression algorithm was used to calculate the half-maximal inhibitory concentration (IC50), and satisfactory prediction accuracy was obtained through 10 times cross-validation (Geeleher et al., 2014). The calculation process was completed by the "pRRophetic" R package.

Construction and Validation of IMRG Prognostic Signature
According to the prognosis-related IMRGs in the univariate Cox regression model, the "glmnet" package in R performed the least absolute shrinkage and selection operator (LASSO) to identify important prognostic IMRGs and select one standard error (SE) above the minimum criteria. The multivariate Cox regression model made it more optimized and practical. Finally, the IMRGscore formula was obtained: Risk score exp Gene1 × coefficient gene1 On the basis of the optimal cut-off of IMRGscore obtained by the "surv_cutpoint" function in R, we divided BLCA patients into high-risk and low-risk groups. With the help of Kaplan-Meier analysis ("survival" package) and receiver operating characteristic (ROC) curve ("timeROC" package), the predictive ability of the prognostic signature was evaluated. The diagnostic accuracy was estimated by the area under the curve (AUC). The same IMRGscore calculation formula, cut-off value and, analysis methods were applied in the GSE13507 cohort to validate the signature.

Establishment and Evaluation of the Nomogram
Nomogram is an intuitive clinical prognosis prediction model integrating a variety of clinicopathological features related to prognosis. We established a nomogram model to provide a more accurate prediction of prognosis for clinical patients based on IMRGscore and clinical pathological characteristics. First, univariate Cox regression analysis was utilized to assess the predicted values of variables. Then further determined the coefficient via multivariate Cox regression analysis. The "rms" R package then established a nomogram for predicting the operating system. Concordance index (C-index) and calibration analysis were applied to estimate the accuracy and consistency. Finally, the clinical application value of the nomogram is evaluated using Decision Curve Analysis (DCA). These analyses were performed with "survival" package.

Statistical Analysis
All statistical analyses were completed with R software (version 4.0.4) in this study. Before establishing and verifying the prognostic signature, the batch differences between the TCGA dataset and GEO datasets were removed through the "sva" package. Wilcoxon rank-sum test was to verify the significance of the difference in two groups. When comparing more than two Frontiers in Cell and Developmental Biology | www.frontiersin.org February 2022 | Volume 10 | Article 810272 groups, the Kruskal Wallis test was selected to verify the difference. Set p-value <.05 as a statistically significant standard.

RESULTS
To describe our research intuitively and systematically, we showed the research process in Figure 1.

Characterization of Iron Metabolism Patterns in BLCA
Through the univariate COX regression analysis (p < .05) of the TCGA-BLCA patients with integrated survival information and cancer tissue expression profile, 95 IMRGs were selected as prognosis-related genes ( Table 1). Then we clustered the TCGA-BLCA cohort by NMF algorithm based on these genes. According to cophenetic coefficients, we decided k = 4 as the best cluster number (Figure 2A). Figure 2B was the NMF matrix heatmap when k = 4, including C1 subtype 89 cases, C2 subtype 141 cases, C3 subtype 91 cases, and C4 subtype 79 cases. Kaplan-Meier survival curves showed that the prognosis of patients in C1 and C3 patterns was better than that of C1 and C2 patterns (p = .020) ( Figure 2C). Figure 2D shows the expression of prognosis-related IMRGs in iron metabolic patterns. We also analyzed the clinicopathological differences among distinct iron metabolism patterns ( Figure 2E). It was found that the proportion of TNM stages, pathologic stage, histologic grade, OS, DSS, and PFI events was dissimilar among patterns, and the incidence of advanced clinicopathological results in C2 and C4 patterns tended to increase.
Through GSVA analysis, we obtained the rich-concentration pathways among iron metabolism patterns ( Figure 2F). We found that C4 pattern was positively related to multiple stromal, carcinogenic, and immune activation related pathways, including TGF-β signaling pathway, WNT signaling pathway, MAPK signaling pathway, JAK-STAT signaling pathway, T cell receptor signaling pathway, chemokine signaling pathway, B cell receptor signaling pathway, cytokine-cytokine receptor interaction, NOD-like receptor signaling pathway, TOLL-like receptor signaling pathway and so on. C2 pattern showed a similar trend to C4 pattern, but C2 was also significantly expressed in a variety of DNA damage repair related pathways. The correlation score of most carcinogenic and immune activation-related pathways was reduced in C1 and C3 patterns. While C3 pattern was also found to exhibit an enriched trend in DNA damage repair related pathways.

Tumor Immune Microenvironment of Iron Metabolism Patterns
In order to investigate whether there are differences in TIME among distinct iron metabolism patterns, we used ESTIMATE and ssGSEA scores for evaluation. ESTIMATE showed that there Frontiers in Cell and Developmental Biology | www.frontiersin.org February 2022 | Volume 10 | Article 810272 were significant differences in the stromal score (p < .001), immune score (p < .001), and ESTIMATE score (p < .001) among the three patterns, of which C4 was the highest, C1 and C3 were the lowest ( Figure 3A). Then we analyzed the infiltration differences in immune cells among iron metabolism patterns. The ssGSEA score suggested that the infiltration of all 22 TIICs in iron metabolism patterns was significantly different, among which, the ssGSEA score of TIICs in C1 and C3 patterns was lower and C2 and C4 patterns were higher ( Figure 3B). And the enrichment trend of immune-related functions in iron metabolism patterns was similar to that of immune cell infiltration ( Figure 3C). Additionally, the expression levels of major histocompatibility complex (MHC) molecules, costimulatory molecules, and adhesion molecules roundly decreased in C1 and C3 patterns ( Figure 3D).

Tumor Somatic Mutation in Iron Metabolism Patterns
The tumorigenesis frequently occurs after the accumulation of gene mutations (Martincorena and Campbell, 2015). It is also reported that tumor mutation burden (TMB) can be used as a potential prognostic indicator for BLCA (Chan et al., 2019). Consequently, we used the "maftools" R package to show the distribution of somatic mutations and the differences of TMB in various iron metabolism patterns. Through the simple nucleotide variation information of TCGA-BLCA, the mutation spectrum and TMB of each sample was obtained. In BLCA samples, the 20 genes with the highest mutation rate were TP53, TTN, KMT2D, MUC16, ARID1A, KDM6A, PIK3CA, SYNE1, RB1, HMCN1, FGFR3, RYR2, KMT2C, MACF1, EP300, FLG, FAT4, STAG2, ATM and OBSCN ( Figures 4A-D). C2 pattern had the highest mutation rate of TP53, while the mutation of TTN and KMT2D mostly happened in C3 pattern. The mutation rates of these three genes in C1 and C4 patterns were significantly reduced. Most gene mutations were missense-mutation. In patients with BLCA, high TMB indicated a better prognosis (hazard ratio [HR] = .65 (.48-.88), p = .005) ( Figure 4E). Additionally, we found that the TMB of C2 and C3 patterns was significantly upper than that of C1 and C4 patterns ( Figure 4F).

Evaluation of Immunotherapy in Iron Metabolism Patterns
Lately, immune checkpoint inhibitors (ICIs) have gradually become the second-line treatment for advanced BLCA. Therefore, we analyzed the expression of some immune checkpoints (PDCD1 (PD-1), CD274 (PD-L1), PDCD1LG2 (PD-L2), LAG3, TIGIT, IDO1, and CTLA4) among different iron metabolism patterns to predict the efficacy of immunotherapy ( Figure 5A). The expression of all immune checkpoints in C2 and C4 patterns was signally higher than that in the other two patterns. This suggested that C2 and C4 might be more suitable for ICIs treatment. However, the high expression level of immune checkpoints may be related to the formation of the immunosuppressive microenvironment (Dunn et al., 2002). This conclusion was confirmed in Figures 5B,C. C4 immune exclusion score was observably higher than other iron metabolism patterns. The immune dysfunction score of C4 pattern also increased signally, while C3 was the lowest. Moreover, we used the TIDE algorithm to evaluate ICIs response, in which the MSI of C1 and the TIDE of C4 were the highest ( Figures 5D,E).

Chemotherapeutics Drugs Response in Iron Metabolism Patterns
Chemotherapeutics drugs are widely used in the treatment of BLCA, including intravesical instillation and systemic chemotherapy. Consequently, we evaluated the IC50 values of  Figures 6A-F). Among the six drugs, C2 patterns showed the lowest IC50 value. In Gemcitabine, Cisplatin, Docetaxel, and Paclitaxel, the IC50 value of C4 pattern was also lower than that of C1 and C3 patterns. As for C3 pattern, the IC50 value in Gemcitabine, Doxorubicin, and Paclitaxel was higher than that of C1 pattern. The above results strongly indicated that C2 was the most sensitive to chemotherapeutic drugs, C4 was the second, while C3 was more insensitive to chemotherapeutic drugs.

Establishment of the IMRG Prognostic Signature in the TCGA-BLCA Cohort
We selected 400 BLCA patients with explicit, non-zero OS and survival status to establish an IMRG signature from the TCGA database. Then used the LASSO Cox regression model to calculate, and selected an SE higher than the minimum standard to further screen the prognostic genes. Finally, through the multivariate COX regression analysis, we obtained the optimal prognostic signature containing 13 IMRGs, including TCIRG1, CTSE, ATP6V0A1, CYP2C8, RNF19A, CYP4Z1, . And according to the optimal cutoff value (cut point = 1.78265), samples were decomposed into low-and high-risk groups. Kaplan-Meier survival analysis showed that the OS of the low-risk group was longer (hazard ratio [HR] = 4.49 (3.29-6.14), p < .001) ( Figure 7C). And the AUCs for the 1-, 3-, and 5-year OS survival rates were .741, .772, and .755, respectively ( Figure 7D). The risk score distribution, survival status graph, and expression profile heatmap were shown in Figures 7E-G. The proportion of patient deaths was observably positively correlated with the IMRGscore. The expression of ATP6V0A1, PLOD1, BMP6, CAST, and SCD were up-regulated in the high-risk group, while TCIRG1, CTSE, CYP2C8, RNF19A, CYP4Z1, YPEL5, IFNG, and ASIC3 were down-regulated.

Confirmation of the IMRG Signature in the GSE13507 Cohort
As the test set, 165 BLCA samples in the GSE13507 cohort were grouped using the same IMRGscore calculation formula and cutoff value of the train set to validate the applicability and stability of the IMRG signature. Consistent with the above conclusion, patients in the low-risk group had a better OS (hazard ratio [HR] = 2.65 (1.52-4.60), p = .001) in the GSE13507 cohort ( Figure 8A). The AUCs for the 1-, 3-, and 5-year OS survival rates were .753, .630, and .552, respectively ( Figure 8B). The conclusions of the risk score distribution, survival status graph, and expression profile heatmap were consistent with the training set ( Figures 8C-E).

Clinical Relevance of the IMRG Signature
To further supplement the clinical application value of the IMRG prognostic signature, we integrated the significant differences in IMRGscore among distinct subgroups of BLCA patients with clinicopathological characteristics ( Figure 9A). Heatmap indicated that the advanced TNM stages, pathologic stage, histologic grade, aging, and worse OS, DSS, and PFI events had an elevated trend in IMRGscore.
Since the significant correlation between signature and clinicopathological stage, we determined whether the IMRGscore was a clinically independent predictor of BLCA patients ( Figures 9B,C). Univariate Cox regression analysis showed that advanced pathologic stage (p < .001), aging (p = .005), and higher IMRGscore (p < .001) were unfavorable factors for OS. After performing the multivariate Cox regression analysis, we confirmed that the IMRGscore was an independent prognostic parameter.

Establishment of a Nomogram Based on the IMRG Signature
According to the above result from univariate Cox regression analysis of the TCGA-BLCA cohort, we established a nomogram model containing pathologic stage, age, and IMRGscore ( Figure 10A). After removing the patients without complete information and the subgroups of variables with few samples, a total of 362 patients were included. We standardized each variable with a score from 0 to 100 and summed the scores to obtain the total number of points for each BLCA patient. The predicted 1-, 3-, and 5-year survival probabilities of each patient were standardized according to the relationship between the positions along with the prognosis and total points axes. The C-index reached 0.694 (95% CI: 0.653-0.735). Figures 10B-D suggested that the nomogram model predicted that the prognosis results of TCGA-BLCA patients would fit well with the actual prognosis results. Besides, DCA curves revealed that the signature provided patients with a stable and significant net benefit in BLCA patients (Figures 10E-G). Then based on the nomogram model, we built the Kaplan-Meier survival curve and the timedependent ROC curves. In the TCGA-BLCA cohort, we divided samples into high-and low-risk groups with the optimal cut-off value (cut point = −.122532). Patients in the high-risk group showed a significantly poor OS (hazard ratio [HR] = 4.22 (2.92-6.10), p < .001, Figure 10H). The AUCs for the 1-, 3-, and 5-year OS survival rates were .764, .769, and .760, respectively ( Figure 10I). Additionally, we verified the nomogram model in the GSE13507 cohort. The cut-off value was consistent with the train set. Kaplan-Meier survival analysis showed that patients in the low-risk group had a better OS than those in the high-risk group (hazard ratio [HR] = 6.47 (2.89-14.49), p < .001,  Figure 10K).

GSEA
To further comprehend the effect of IMRGs expression on the biological characteristics of BLCA, we carried on GSEA analysis in IMRGscore-defined groups ( Figures 11A,B). The Kyoto Encyclopedia of Genes and Genomes (KEGG) results revealed that in the high-risk group, the main enrichment pathways were ECM receptor interaction, regulation of actin cytoskeleton, MAPK signaling pathway, WNT signaling pathway, pathways in cancer. While the low-risk group is mainly concentrated in allograft rejection, asthma, primary immunodeficiency, and so on. Furthermore, Figures 11C,D showed the enrichment of the high-and low-risk groups in the Gene Ontology biological process (GOBP). We found that the low-risk group was enriched in multiple immune functions, such as activation of immune response, adaptive immune response, B cell-mediated immunity, and so on.

Correlation Between Tumor Immune Microenvironment and IMRG Signature
As we mentioned before, iron metabolism is closely related to TIME. Based on the CIBERSORT algorithm, we calculated the proportion of 22 TIICs in each TCGA-BLCA sample ( Figure 12A). After selecting samples with significant immune cell fraction results (p < .05), 195 samples were included in the difference analysis, including 141 cases in the low-risk group and 54 cases in the high-risk group. Then the difference in the proportion of TIICs between the IMRGscore-defined groups was shown in Figure 12B. It was found that the fraction of CD8 T cells, activated CD4 memory T cells, follicular helper T cells (TFH) and regulatory T cells (Treg) in the low-risk group was significantly higher. In contrast, the fraction of M0 macrophages was upper in the high-risk group. Furthermore, high levels of CD8 T cells (p = .004), activated CD4 memory T cells (p = .013) and TFH (p = .041) were observably associated with better OS (Figures 12C-E), while increased M0 macrophages (p = .035) indicated a poor OS ( Figure 12F).

DISCUSSION
Bladder cancer is a heterogeneous malignancy. Patients with BLCA generally show different prognoses because of the molecular discrepancy (Knowles and Hurst, 2015;Guo and Czerniak, 2019). At present, it is generally believed that clinical or pathological stages were insufficient to predict the prognosis of patients with BLCA (Konety, 2006;Rosenberg et al., 2013). Therefore, it is essential to develop a more accurate and efficient model to predict the survival prognosis for patients. In recent years, some studies have found iron involvement in the appearance and progression of cancers. Abundant researches have revealed that iron metabolism is involved in the entire process of cancer progression. Murata M. demonstrated that iron released from the damaged transferrin could mediate the Fenton reaction and produce ROS, which contributes to the carcinogenic process in multiple ways (Murata, 2018). A crosssectional study found that the serum iron concentration in patients with BLCA was lower than that in the control group (Mazdak et al., 2010). However, studies on the potential function of iron metabolism in the treatment and prognosis of BLCA are scarce.
In this study, patients in the TCGA-BLCA cohort were divided into four iron metabolism patterns based on the expression of prognosis-related IMRGs. Survival analysis suggested that C1 and C3 had a better prognosis. Our results also showed that C2 and C4 patterns have higher enrichment scores in multiple carcinogenic and immune activation-related pathways. For instance, abnormal activation of NOD-like receptors (NLRs) occurs in various cancers, coordinates the tumor immune microenvironment, and promotes angiogenesis, cancer cell stem cells, and chemotherapy resistance, thereby enhancing tumor risk . Toll-like receptor (TLR) is a transmembrane pattern recognition receptor that detects and defends microbial pathogens through the innate immune response (Brennan and Gilmore, 2018). The activation of the Notch signaling pathway can be seen in most components of the tumor microenvironment (TME), such as angiogenesis, tumor stem cell maintenance, immune infiltration, or therapeutic resistance (Meurette and Mehlen, 2018). Besides, C2 and C4 iron metabolism patterns were highly infiltrated by immune cells, Frontiers in Cell and Developmental Biology | www.frontiersin.org February 2022 | Volume 10 | Article 810272 13 and the expression of MHC genes was highly elevated. These characteristics are consistent with an immune-inflamed phenotype. On the other hand, C1 and C3 are suspected to be immune-desert phenotypes due to lack of immune infiltration and antigen presentation. Additionally, we found that multiple immune checkpoints (PD-1, PD-L1, PD-L2, LAG3, TIGIT, IDO1, and CTLA4) were highly expressed in C2 and C4 patterns, which might indicate that immunotherapy has a better effect on them. Furthermore, studies have shown that the high-level expression of immune checkpoint genes might be a stimulative to the immunosuppressive microenvironment and led to the immune escape of tumor cells (Dunn et al., 2002). The TIDE algorithm confirmed this conclusion. Therefore, we speculate that the reason why C2 and C4 patterns do not show a better prognosis of immunoinflammatory phenotype is that the antitumor effect based on the activated immune pathway and high infiltration level of T cells were eliminated by the formation of the immunosuppressive microenvironment. The above evidence proved that iron metabolism is of great significance in shaping various TME landscapes in BLCA.
Somatic mutation is not only the driving factor of cancer, but also the guiding basis for diagnosis and treatment. The three genes with the highest mutation rate in BLCA were TP53, TTN, and KMT2D. The mutation rate of TP53 in C2 iron metabolism pattern was the highest, while in C1 pattern was the lowest. In C2 and C3 patterns, the incidence of TTN and KMT2D mutations was significantly higher. Detection of TP53 mutation was conducive to estimating the high risk of early lesions (Olivier et al., 2010). Single TTN gene mutation often indicated high TMB (Oh et al., 2020). The mutation of epigenetic regulator KMT2D was a biomarker of poor prognosis in some cancers (Ferrero et al., 2020). Additionally, we found that high TMB suggested that patients with BLCA had a better prognosis through the TCGA database. Consistently, the TMB of C1 and C4 patterns was observably decreased. These results show that iron metabolism has a complex interaction with somatic mutation. Chemotherapy is still one of the main treatments for BLCA. According to the latest guidelines of the European Association of Urology (EAU), all muscle-invasive bladder cancer (MIBC) patients with physical conditions can apply platinum-based neoadjuvant chemotherapy before operation (Witjes et al., 2021). This study investigated the efficacy of six common chemotherapeutic drugs on iron metabolism mode, including Gemcitabine, Cisplatin, Docetaxel, Mitomycin-C, Doxorubicin, and Paclitaxel. The results showed that C2 pattern was the most sensitive to these chemotherapeutic drugs, while C3 pattern was the most insensitive. This provided a reference basis for the selection of clinical chemotherapy drugs.
Due to the effect of iron metabolism on the tumorigenesis and progression of cancer, it is of great importance to establish an IMRG signature for predicting the prognosis of patients with BLCA. Herein, we applied stepwise regression analysis to compose a clinical prognostic signature for BLCA patients with 13 IMRGs (TCIRG1, CTSE, ATP6V0A1, CYP2C8, RNF19A, CYP4Z1, YPEL5, PLOD1, BMP6, CAST, SCD, IFNG, and ASIC3). A test set was utilized to confirm its accuracy and stability. Moreover, the IMRGscore was elevated in multiple advanced clinicopathological stages. And multivariate Cox regression analysis verified that IMRGscore was an independent prognostic index of BLCA patients. Finally, we combined some clinicopathological features, including pathologic stage, age, and IMRGscore to construct a nomogram that accurately predicted the prognosis of patients with BLCA. The accuracy and clinical contributions were verified by calibration analysis and DCA.
Several studies have shown that these 13 IMRGs are closely related to cancers, and some of these genes have been confirmed about BLCA. TCIRG1, one of the V-ATPase subunits, is abnormally overexpressed in patients with recurrent hepatocellular carcinoma, and enhances the ability of metastasis by regulating the growth, death, and epithelial to mesenchymal transition of cancer cells (Yang et al., 2018). A retrospective study suggests that CTSE can be used as an independent prognostic marker for NMIBC, so as to guide the treatment of patients (Lin et al., 2001). ATP6V0A1was demonstrated that it could enhance the fusion of autophagosomes and lysosomes, up-regulate autophagy volume accumulation, and finally induce autophagic cell death (Hsin et al., 2012). CYP2C8 can be regulated by GAS5/miR-382-3p in hepatocellular carcinoma and play an anticancer role (Li and Chen, 2020). RNF19A was confirmed to be overexpressed in non-small cell lung cancer, which plays a carcinogenic role by destroying the function of p53 (Cheng et al., 2021). CYP4Z1 was confirmed to be highly expressed in BLCA and positively correlated with the progression of histologic grade and pathologic stage (Al-Saraireh et al., 2021). YPEL5 was found to be inhibited by METTL3-m6A (N6-methyladenosine)-YTHDF2 axis in colorectal cancer, promoting the growth and metastasis of tumor (Zhou et al., 2021). The expression of PLOD1 was an independent prognostic factor in BLCA patients, and downregulated by inhibitor could significantly reduce the invasiveness of BLCA cells (Yamada et al., 2019). The expression of BMP6, a key endogenous regulator of iron metabolism, was affected by Med19, which could promote bone metastasis and invasiveness of bladder cancer (Andriopoulos Jr et al., 2009;Wen et al., 2013). Public clinical data also confirm that BMP6 is a prognostic marker for bladder cancer (Yuen et al., 2012). Calpastatin (CAST) is involved in many important physiological processes, including cell cycle, ECM, cancer cell proliferation, metastasis, and apoptosis (Nian and Ma, 2021). SCD can protect cancer cells from oxidative stress and ferroptosis through mediated lipogenesis in prostate cancer with over-activation of PI3K-AKT-mTOR signaling (Yi et al., 2020). SCD has also been shown to reduce proliferation and invasion of BLCA cells when inhibited (Piao et al., 2019). Iron metabolism can affect innate immune response by affecting IFNG mediated immune response pathway in macrophages (Nairz et al., 2014). And IFNG was demonstrated to inhibit the activity of bladder cancer stem cells (Qiu et al., 2020). ASIC3 has an H+ gating function, which promotes the acid-induced epithelial-mesenchymal transition in pancreatic cancer cells (Zhu et al., 2017). GSEA analysis showed that several cancer-related pathways were enriched in the high-risk group. The unique biochemical and biophysical properties of ECM when it is dysregulated are the key drivers of cancer progression (Walker et al., 2018). The MAPK signaling pathway is considered to be related to cell proliferation, differentiation, migration, aging, and apoptosis (Sun et al., 2015). The Wnt signaling pathway is an important driving factor to maintain tissue development and homeostasis. Abnormal Wnt signaling will cause the occurrence and progression of many cancers by affecting cancer stem cells (Duchartre et al., 2016). Furthermore, CIBERSORT analysis indicated that several TIICs with differential infiltration levels had significant correlations between the prognosis of BLCA patients, and played a regulatory role in the balance of iron metabolism. After being activated by immunotherapy, CD8 T cells can enhance ferroptosis-specific lipid peroxidation and increase ferroptosis in tumor cells, which contributes to the antitumor effect . Macrophages are the regulatory hub of iron metabolism. Macrophages phagocytize and degrade aging and damaged erythrocytes to recover iron, and also have the ability to release iron. The accumulation of M2 macrophages is often associated with poor prognosis, and M2 macrophages possess iron release characteristics (Recalcati et al., 2010). Therefore, the accumulation of M2 macrophages may further aggravate the disorder of iron metabolism.
Our research has obvious advantages in the study of iron metabolism characteristics and the clinical application of BLCA. First, we divided BLCA patients into distinct iron metabolism landscapes to further confirm the relationship between iron metabolism and TME and somatic mutations. At the same time, we also verified that iron metabolism landscapes had guiding significance for chemotherapy drugs and immunotherapy selection. Next, we established the IMRG prognostic signature and proved IMRGscore was an independent prognostic factor for BLCA patients, and it also had the predictive ability for clinicopathological characteristics. At present, our research is still insufficient. First, it is a retrospective study. Deviations in variables such as clinicopathological characteristics of patients most likely exist. Second, our prognostic signature depends on gene expression and does not take into account the effects of gene mutation, methylation, or other factors on the prognosis of BLCA patients. Finally, the prognostic signature can be incorporated into large sample prospective studies to further verify its clinical value.

DATA AVAILABILITY STATEMENT
Raw data for this study were generated at the TCGA database with the cancer type of BLCA. The datasets used and/or analyzed during the current study are available from the GEO database (GSE13507). Derived data supporting the findings are available from the corresponding author (WS) on reasonable request.

AUTHOR CONTRIBUTIONS
XL and WS designed and guided the work. XS and SX participated in data collecting, data processing, figures preparation, and manuscript writing. YZ, JM, and CD contributed to statistical analysis. KC, LC, and FL contributed to manuscript writing and article publishment. ZL and TW contributed to manuscript draft writing. JL revised the manuscript critically. All authors provided critical advice for the final manuscript.