A Novel Ferroptosis-related Gene Signature for Overall Survival Prediction in Patients with Bladder cancer

Background Bladder cancer (BC) is a highly heterogeneous disease, which makes the prognostic prediction challenging. Ferroptosis is related to a variety of biological pathways, including those involved in the metabolism of amino acids, lipids, and iron. However, the prognostic value of ferroptosis-related genes in BC remains to be further elucidated. Methods In this study, the mRNA expression proles and corresponding clinical data of BC patients were downloaded from public databases. The least absolute shrinkage and selection operator (LASSO) Cox regression model was utilized to construct a multigene signature and validated it. Results Our results showed 12 differentially expressed genes (DEGs) were correlated with overall survival (OS) in the univariate Cox regression analysis (all adjusted P< 0.05). A 9-gene signature was constructed to stratify patients into two risk groups. Patients in the high-risk group showed signicantly reduced OS compared with patients in the low-risk group (P < 0.001). The risk score was an independent predictor for OS in multivariate Cox regression analyses (HR> 1, P< 0.01). Receiver operating characteristic (ROC) curve analysis conrmed the signature's predictive capacity. Functional analysis revealed that immune-related pathways were enriched, and immune status were different between two risk groups, especially in humoral immune response process. Conclusion In conclusion, a novel gene Comparison of the scores between different risk groups in the TCGA cohort. The scores of 16 immune cells (a) and 13 immune-related functions (b) are displayed in boxplots. CCR, cytokine-cytokine receptor. Adjusted P values were showed as: ns, not signicant; *, P< 0.05; **, P< 0.01;


Introduction
Bladder cancer (BC) is a common cancer with 81,400 new cases in 2020 [1],and death rates remained stable between 2007 and 2016. Early stages of disease are often treated with transurethral resection, and/or complete cystectomies, and/or neoadjuvant or adjuvant cisplatin-based chemotherapy. Treatment selection is often decided by tumor stage, risk factors, performance status and comorbidities. However approximately 5% of patients present with metastatic BC at initial diagnosis, a large portion relapses or progresses to advanced stages following treatment for localized disease with a 5-year relative survival of 4.6%. Because of no speci c tumor markers of BC makes prognostic prediction challenging, and it's a very easy recurrence. There is a necessary to develop of novel prognostic models and nd new therapeutic target.
Ferroptosis is related to a variety of biological pathways, including those involved in the metabolism of amino acids, lipids, and iron [2,3]. At present, ferroptosis is known to be closely related to the occurrence, development, and treatment of cancer and to play a dual role in tumor immunity. Ferroptosis promotes tumor cell death to facilitate anti-tumor immunity; on the other side, it helps neighboring tumor cells escape immune-mediated killing and promotes tumor growth [4,5]. Previous studies have reported that ferroptosis plays a vital role, such as CISD1 [6] and TP53 gene [7]. On the other hand, other ferroptosisrelated genes, such as Rb[8], NRF2 [9] might protect tumor to avoid ferroptosis. But there were rarely papers about ferroptosis in bladder cancer. As we all know, Martin-Sanchez et al. reported intracellular iron concentration relation with bladder tumor proliferation, found that ion concentration to drop with bladder cancer cell proliferation in 1993 [10]. However, the writers didn't make further research, and whether these ferroptosis-related genes are correlated with BC patient prognosis remains unknown. In the present study, we rstly to research whether these ferroptosis-related genes are correlated with BC, through mRNA expression and corresponding clinical data from public databases. And we constructed a prognostic multigene signature using ferroptosis-related differentially expressed genes (DEGs) and validated it.

Data collection
The RNA sequencing and clinical data of 433 BC patients were downloaded from the TCGA website (https://portal.gdc.cancer.gov/repository). The gene expression pro les were normalized using the scale method provided in the "limma" R package. Normalized read count values were used. The current research follows the TCGA data access policies and publication guidelines.
Construction and validation of a prognostic ferroptosis-related gene signature The "limma" R package was used to identify the differentially expressed genes (DEGs) between tumor tissues and adjacent nontumorous tissues (|log2FC| ≥ 1, FDR < 0.05). Univariate Cox analysis of overall survival (OS) was performed to screen ferroptosis-related genes with prognostic values. P values were adjusted by Benjamini & Hochberg (BH) correction. Use the R package "Venn" screen intersection of iron death related genes with differential and prognostic value. An interaction network for the overlapping prognostic DEGs was generated by the STRING database (version 11.0) [11]. And then all samples' data were divided into train group and test group averagely. The LASSO-penalized Cox regression analysis was applied to construct a prognostic model using train group to minimize the risk of over tting [12,13].
The LASSO algorithm was used for variable selection and shrinkage with the "glmnet" R package. The independent variable in the regression was the normalized expression matrix of candidate prognostic DEGs, and the response variables were overall survival and status of patients in the train group. Penalty parameter (λ) for the model was determined by tenfold cross-validation following the minimum criteria (i.e. the value of λ corresponding to the lowest partial likelihood deviance). The risk scores of the patients were calculated according to the normalized expression level of each gene and its corresponding regression coe cients. The formula was established as follows: score= e sum (each gene' s expression × corresponding coe cient) . The patients were strati ed into high-risk and low-risk groups based on the median value of the train risk score. Based on the expression of genes in the signature, PCA was carried out with the "prcomp" function of the "stats" R package. Besides, t-SNE were performed to explore the distribution of different groups using the "Rtsne" R package. For the survival analysis of each gene, the optimal cut off expression value was determined by the "surv_cutpoint" function of the "survminer" R package. The "survivalROC" R package was used to conduct time-dependent ROC curve analyses to evaluate the predictive power of the gene signature. And then validated all in the test group.

Functional enrichment and immune status analysis
The "clusterPro ler" R package was utilized to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on the DEGs (|log2FC| ≥ 1, FDR < 0.05) between the highrisk and low-risk groups. P values were adjusted with the BH method. The in ltrating score of 16 immune cells and the activity of 13 immune-related pathways were calculated with single-sample gene set enrichment analysis (ssGSEA) [14] in the "gsva" R package. The annotated gene set le is provided in Supplementary Table S2.

Statistical analysis
Student's t-test was used to compare gene expression between tumor tissues and normal tissues.
Differences in proportions were compared by the Chi-squared test. Mann-Whitney test with P values adjusted by the BH method was used to compare the ssGSEA scores of immune cells or pathways between the high risk and low risk group. The OS between different groups was compared by Kaplan-Meier analysis with the log-rank test. Univariate and multivariate Cox regression analyses were implemented to identify independent predictors of OS. All statistical analyses were performed with R software (Version 3.6.1). If not speci ed above, a P value less than 0.05 was considered statistically signi cant.

Results
The ow chart of this study is shown in Fig. 1. A total of 433 BC patients from the TCGA-LIHC cohort were nally enrolled. The detailed clinical characteristics of these patients are summarized in Table 1.
Identi cation of prognostic ferroptosis-related DEGs in the TCGA cohort 61 ferroptosis-related genes were differentially expressed between tumor tissues and adjacent nontumorous tissues (Supplementary Table S3), and 12 of them were correlated with OS in the univariate Cox regression analysis (Fig. 2a). A total of 12 prognostic ferroptosis-related DEGs were preserved ( Fig.2b-c). The interaction network among these genes indicated that SLC2A12 and SLC3A2 were the hub genes (Fig.2d). The correlation between these genes is presented in Fig. 2e.

Construction of a prognostic model in the train group
All samples data were divided into train group (n=200) and test group (n=200) averagely (Supplementary   Table S5 and S6). LASSO Cox regression analysis was applied to establish a prognostic model using the expression pro le of the 12 genes mentioned above in train group. After LASSO Cox regression analysis, the optimal 9 genes and their corresponding correlation coe cients were obtained (Supplementary Table  S4). Survival analyses, according to the optimal cut-off expression value of each gene, indicated that high expression of these genes all correlated with a poor prognosis (all adjusted P<0.05, Fig. S2). The risk score was calculated as follows: e (0.001 * expression level of CAPG + 0.059 * expression level of CDO1 + 0.609 * expression level of DRD5 + 0.001 * expression level of SCD + 0.003 * expression level of SLC3A2 + 0.003 * expression level of TFRC -0.072 * expression level of SLC38A1 -0.047 * expression level of TP63 -0.370 * expression level of ZNF419). The patients were strati ed into a high-risk group (n=100) or a low-risk group (n=100) according to the median cut-off value (Fig. 3a). As shown in Fig. 3b patients with high risk had a higher probability of death earlier than those with low risk. PCA and t-SNE analysis indicated the patients in different risk groups were distributed in two directions ( Fig. 3c-d). Consistently, the Kaplan-Meier curve indicated that patients in the high-risk group had a signi cantly worse OS than their low-risk counterparts (Fig. 3e, P< 0.001). The predictive performance of the risk score for OS was evaluated by time-dependent ROC curves, and the area under the curve (AUC) reached 0.712 at 1 year, 0.661 at 2 years, and 0.689 at 3 years (Fig. 3f).

Validation of prognostic model in test group
Using test group to verify the robustness of the prognostic model constructed. The same risk scoring formula was used to calculate, and the calculated median value (train group) was divided into high-risk group and low-risk group, as shown in Figure 4a. Similarly, patients in the high-risk group had a higher rate of early death than patients in the low-risk group, as shown in Figure 4b. PCA and t-SNE analysis indicated the patients in different risk groups were distributed in two directions (Fig. 4c-d). Consistently, the Kaplan-Meier curve indicated that patients in the high-risk group had a signi cantly worse OS than their low-risk counterparts (Fig. 4e, P< 0.01). The predictive performance of the risk score for OS was evaluated by time-dependent ROC curves, and the area under the curve (AUC) reached 0.682 at 1 year, 0.693 at 2 years, and 0.706 at 3 years (Fig. 4f).
Independent prognostic value of the 9-gene signature Univariate and multivariate Cox regression analyses were carried out among the available variables to determine whether the risk score was an independent prognostic predictor for OS. In univariate Cox regression analyses, the risk score was signi cantly associated with OS (HR= 3.599, 95% CI =2.383-5.437, P< 0.001). After correction for other confounding factors, the risk score still proved to be an independent predictor for OS in the multivariate Cox regression analysis (HR =3.011, 95% CI = 1.956-4.635, P<0.001) ( Fig. 5a, b).

Functional analyses in the TCGA cohort
To elucidate the biological functions and pathways that were associated with the risk score, the DEGs between the high-risk and low-risk groups were used to perform GO enrichment and KEGG pathway analyses. As expected, DEGs were enriched in several ferroptosis-related molecular functions(P. adjust < 0.05, Fig.6a, c). Interestingly, the DEGs from the TCGA cohort were also obviously enriched in many immune-related biological processes, including humoral immune response, complement activation, classical pathway, humoral immune response mediated by circulating immunoglobulin, cell recognition, phagocytosis, recognition, immunoglobulin complex, circulating, antigen binding, immunoglobulin receptor binding, cytokine activity, chemokine activity, and chemokine receptor binding (P. adjust < 0.05, Fig. 6c). KEGG pathway analyses also indicated that the Viral protein interaction with cytokine and cytokine receptor, cytokine-cytokine receptor interaction pathway, Chemokine signaling pathway, Complement and coagulation cascades, Intestinal immune network for IgA production and Natural killer cell mediated cytotoxicity was enriched in both cohorts (P < 0.05, Fig. 6b, d).
To further explore the correlation between the risk score and immune status, we quanti ed the enrichment scores of diverse immune cell subpopulations, related functions, or pathways with ssGSEA. Interestingly, the score of most immune cells were signi cantly different between the low risk and high risk group in the TCGA cohort, except NK cells, T helper cells and Th 2 cells (all adjusted P< 0.05, Fig. 7a). Moreover, the score of most immune function were signi cantly different between the low risk and high risk group.
lower in the high risk group, while the type II IFN response were no signi cantly (adjusted P< 0.05, Fig. 7ab).

Discussion
Ferroptosis is a form of necrotic cell death marked by the oxidative modi cation of phospholipid membranes via an iron-dependent mechanism [15]. Although a few previous ferroptosis related studies in bladder cancer[16], the research on ferroptosis is in its initial stage. And the mechanism of ferroptosis in the development, treatment and prognosis of urinary tumors has not been fully elaborated. In the current study, we systematically investigated the expression of 259 ferroptosis-related genes in BC tumor tissues and their associations with OS. A novel prognostic model integrating 9 ferroptosis-related genes was rstly constructed and validated it. Functional analyses revealed that immune-related pathways were enriched.
The prognostic model proposed in the present study was composed of 10 ferroptosis-related genes (CAPG, CDO1, DRD5, SCD, SLC3A2, TFRC, SLC38A1, TP63, ZNF419). These genes could be roughly classi ed into three categories, including iron metabolism (TFRC), lipid metabolism (SCD), glutathione metabolism (CAPG, CDO1, DRD5, SLC3A2, SLC38A1, TP63). CAPG (Capping actin protein, gelsolin like) encodes a member of the gelsolin/villin family of actin-regulatory proteins. May play a role in regulating cytoplasmic and/or nuclear structures through potential interactions with actin. It has been reported CAPG was a ferroptosis down-regulated factors which stimulate and suppress ferroptosis by affecting the synthesis of GSH [17]. CDO1(cysteine dioxygenase type 1) also a ferroptosis down-regulated factor. Increasing CDO1 activity can prevent cytotoxicity from elevated cysteine levels. However, cellular cysteine is an indispensable substrate for the synthesis of GSH [18]. So CDO1 competitively limits the availability of cysteine to synthesize GSH and in uences the process of ferroptosis [19]. DRD5(dopamine receptor D5) a subtype of the dopamine receptor. Dopamine can reduced erastin-induced ferrous iron accumulation, glutathione depletion, and malondialdehyde production [20]. Moreover, dopamine promoted dopamine receptor D5 gene expression. DRD5 can protect cell to some extent. SCD (stearoyl-CoA desaturase) encodes an enzyme involved in fatty acid biosynthesis, primarily the synthesis of oleic acid. The protein belongs to the fatty acid desaturase family and plays an important role in regulating the expression of genes that are involved in lipogenesis and in regulating mitochondrial fatty acid oxidation. There was study that have reported inhibition of SCD1 induces ferroptotic cell death [21]. SLC3A2(solute carrier family 3 member 2) is a member of the solute carrier family and encodes a cell surface, transmembrane protein of glutamate-cystine antiporter system xc- [22]. TFRC (Transferrin receptor protein 1) uptake of iron occurs via receptor-mediated endocytosis of ligand-occupied transferrin receptor into specialized endosomes. Endosomal acidi cation leads to iron release. RNAi of transferrin receptor (TfR) inhibited ferroptosis. TP63(Tumor protein 63) acts as a sequence speci c DNA binding transcriptional activator or repressor. May be required in conjunction with TP73/p73 for initiation of p53/TP53 dependent apoptosis in response to genotoxic insults and the presence of activated oncogenes. Overexpression protects cells from ferroptosis-inducing agents [23]. ZNF419(Zinc nger protein 419) has been reported up-regulated (>= 2 fold) in erastin-treated samples and it may promote ferroptosis. Relevant experimental con rmation is needed.
Although the mechanisms underlying tumor susceptibility to ferroptosis have been an intense area of research in the past few years, the potential modulation between tumor immunity and ferroptosis remains elusive. Based on the DEGs between different risk groups, we performed GO analyses and unexpectedly discovered that many immune-related biological processes and pathways were enriched. It is reasonable to assume that ferroptosis may have a close connection with tumor immunity. A typical hallmark of tumors is immune dysregulation, and the interaction between ferroptosis and lipid metabolism is critical in tumor immunomodulation. ferroptosis cells affect immune cells by releasing lipid metabolites, and the mobilization of immune cells to tumors is controlled by various chemokines, cytokines, and soluble metabolites (such as lipids and nucleic acids). However, it is not yet clear whether ferroptosis cells can produce this signal, but some studies have shown that macrophages can clear the cells undergoing ferroptosis, supporting the existence of the signal [24]. Immune cells also regulate ferroptosis in tumor cells. During tumor invasion, initial CD8 + T cells differentiate into effector CD8 + T cells, and are further activated into cytotoxic CD8 + T cells and memory CD8 + T cells, which play a local targeting role in the tumor. Tumor immunotherapy can enhance the role of CD8 + T cells in the tumor microenvironment.
Immunotreaty-activated CD8 + T cells down-regulated the expression of SLC3A2 and SLC7A11, two subunits of System XC -, through interferon gamma (IFN-γ) and other cytokines, which resulted in the failure of GSH in tumor cells to remove lipid peroxides through GPX4, and promoted lipid peroxidation and ferroptosis of tumor cells [22]. ferroptosis directly affects the phenotype and function of immune cells. For example, iron-dead tumor cells release the cancer-causing KRAS protein, which can be swollen tumor -associated macrophages phagocyte and drive their polarization [25]. Interestingly, the contents of the humoral immune response process were signi cantly different between the low risk and high risk group in this study.
There are several limitations of this study. First, our prognostic model was both constructed and validated with retrospective data from public databases. More prospective real-world data are warranted to verify its clinical utility. Second, the intrinsic weakness of merely considering a single hallmark to build a prognostic model was inevitable, as many prominent prognostic genes in BC might have been excluded. In addition, it should be emphasized that the links between the risk score and immune activity have not yet been experimentally addressed.

Conclusion
In summary, this study de ned a novel prognostic model of 9 ferroptosis-related genes for BC. This model proved to be independently associated with OS in both the derivation and validation cohorts, providing insight into the prediction prognosis. However the underlying mechanisms between ferroptosis-related genes and tumor immunity in BC remain poorly understood and warrant further investigation.