Construction and validation of a bladder cancer risk model based on autophagy-related genes

Autophagy has an important association with tumorigenesis, progression, and prognosis. However, the mechanism of autophagy-regulated genes on the risk prognosis of bladder cancer (BC) patients has not been fully elucidated yet. In this study, we created a prognostic model of BC risk based on autophagy-related genes, which further illustrates the value of genes associated with autophagy in the treatment of BC. We first downloaded human autophagy-associated genes and BC datasets from Human Autophagy Database and The Cancer Genome Atlas (TCGA) database, and finally obtained differential prognosis-associated genes for autophagy by univariate regression analysis and differential analysis of cancer versus normal tissues. Subsequently, we downloaded two datasets from Gene Expression Omnibus (GEO), GSE31684 and GSE15307, to expand the total number of samples. Based on these genes, we distinguished the molecular subtypes (C1, C2) and gene classes (A, B) of BC by consistent clustering analysis. Using the genes merged from TCGA and the two GEO datasets, we conducted least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression analysis to obtain risk genes and construct autophagy-related risk prediction models. The accuracy of this risk prediction model was assessed by receiver operating characteristic (ROC) and calibration curves, and then nomograms were constructed to predict the survival of bladder cancer patients at 1, 3, and 5 years, respectively. According to the median value of the risk score, we divided BC samples into the high- and low-risk groups. Kaplan–Meier (K-M) survival analysis was performed to compare survival differences between subgroups. Then, we used single sample gene set enrichment analysis (ssGSEA) for immune cell infiltration abundance, immune checkpoint genes, immunotherapy response, gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis, and tumor mutation burden (TMB) analysis for different subgroups. We also applied quantitative real-time polymerase chain reaction (PCR) and immunohistochemistry (IHC) techniques to verify the expression of these six genes in the model. Finally, we chose the IMvigor210 dataset for external validation. Six risk genes associated with autophagy (SPOCD1, FKBP10, NAT8B, LDLR, STMN3, and ANXA2) were finally screened by LASSO regression algorithm and multivariate Cox regression analysis. ROC and calibration curves showed that the model established was accurate and reliable. Univariate and multivariate regression analyses were used to verify that the risk model was an independent predictor. K-M survival analysis indicated that patients in the high-risk group had significantly worse overall survival than those in the low-risk group. Analysis by algorithms such as correlation analysis, gene set variation analysis (GSVA), and ssGSEA showed that differences in immune microenvironment, enrichment of multiple biologically active pathways, TMB, immune checkpoint genes, and human leukocyte antigens (HLAs) were observed in the different risk groups. Then, we constructed nomograms that predicted the 1-, 3-, and 5-year survival rates of different BC patients. In addition, we screened nine sensitive chemotherapeutic drugs using the correlation between the obtained expression status of risk genes and drug sensitivity results. Finally, the external dataset IMvigor210 verified that the model is reliable and efficient. We established an autophagy-related risk prognostic model that is accurate and reliable, which lays the foundation for future personalized treatment of bladder cancer.


Introduction
BC is the most common tumor in the urinary system. It is reported that BC has a high recurrence rate, with approximately half of patients recurring and developing metastases after radical surgery (Yan et al. 2021;Wu et al. 2020). The existing known histological grading and the depth of tumor invasion have a great impact on its prognosis. Depending on the depth of tumor invasion, BC is divided into non-muscle invasive bladder cancer (NMIBC) and muscle invasive bladder cancer (MIBC), of which 70% is NMIBC and 30% is MIBC (Yang et al. 2021). Traditional surgery combined with chemoradiotherapy still remains highly insufficient for patients with metastatic tumors, and the sensitivity of treatment is variable. Therefore, the discovery of biomarkers for early diagnosis and prediction of prognosis is essential to improve the outcome of BC.
Autophagy is a fundamental evolutionary catabolic process in cells, which maintains cellular homeostasis by degrading and recycling damaged organelles through autophagosomes. However, autophagy is a double-edged sword, and the autophagic process in various conditions of stress also plays a crucial role in the development and progression of cancer, as it not only inhibits tumorigenesis but also facilitates the proliferation of cancer cells, thus promoting tumor development .
TNM staging of tumors is now widely used to predict the prognosis of patients with BC. However, these traditional treatment strategies are not ideal and have a high rate of metastasis and recurrence. Bacillus Calmette-Guerin (BCG) vaccine is a common immunochemotherapy drug used for patients who have a high risk of recurrence of NMIBC, but challenges still exist (Tran et al. 2021). Currently, immunotherapy has made some progress in cancer treatment, and inhibitors of immune checkpoint genes are increasingly used, and new evidence suggests that checkpoint inhibitors programmed cell death 1 ligand 1 (PD-L1), programmed cell death 1 (PD-1), and cytotoxic T-lymphocyte antigen 4 (CTLA-4) show promising results in some cancers (Ma et al. 2020). According to the research (Folkerts et al. 2019;Schaaf et al. 2019), the immune microenvironment plays a key role in tumor progression, immune escape, and immunotherapeutic response. Therefore, a better understanding of the immune microenvironment may help improve BC treatment strategies. TMB is also a novel feature of cancer. As our knowledge of the molecular biology of BC improves, more mechanisms that promote bladder development are being discovered, which allows us to explore new predictive models and biomarkers that might be used to develop personalized treatment plans for BC patients. In this study, our multi-omics model constructed by exploring the role of autophagy-related genes, tumor microenvironment (TME), TMB, and HLA in BC may help provide a deeper understanding of the mechanisms of BC immunotherapy and offer new therapeutic targets for BC.

Data acquisition and processing
We downloaded autophagy-related genes from the Human Autophagy Gene Database (http:// www. autop hagy. lu/ index. html) and samples from the TCGA database for BC, respectively. Using differential analysis of cancer versus neighboring cancers and univariate Cox regression analysis, setting false discovery rate (FDR) < 0.05 and |LogFC|> 1.5 as thresholds for screening, we obtained 44 autophagyassociated genes differentially expressed in tumor tissue and surrounding tumors. The datasets, GSE31684 and GSE13507, were downloaded from the GEO database (http:// www. ncbi. nlm. nih. gov/ geo/), and these two datasets were merged with the TCGA database using the "ComBat" function, and the "sva" package was used to eliminate nonbiologically biased batch effects, and finally filtered out 39 differential autophagy-related genes.

Identification of BC subclasses and signature construction
We conducted a consistent clustering analysis using the cancer subtypes package based on the 39 autophagyrelated genes that had been obtained, and obtained the molecular subtypes C1 and C2 for BC. Then, using the TCGA and GSE31684/GSE15307 databases, we carried out the LASSO and multivariate Cox regression analysis to obtain six risk prediction genes (SPOCD1, FKBP10, NAT8B, LDLR, STMN3, and ANXA2) and constructed a risk prediction model with the following formula: Risk score = nΣi = Exp(mRNA 1)*coefficient 1 + Exp(miRNA 2)*coefficient 2 + … + Exp(miRNA n)*coefficient n.
Next, we randomly distributed the samples merged from the TCGA and two datasets of GEO into the training and test sets in the ratio of 7:3. A comparison of the training and test datasets is shown in Table 1 and then performed the internal validation of the models separately. According to the median risk score, we distinguished BC patients into the high-risk and low-risk groups. We used the K-M survival function to assess the difference in overall survival across risk groups. Subsequently, we explored genes within survival differences and obtained gene groupings by cluster analyses A and B. Finally, ROC and calibration curves were used to assess the accuracy of the model. Meanwhile, univariate and multivariate Cox regression analyses showed that risk score was an independent predictor.

Pathway and function enrichment analysis
To investigate the potential mechanisms that affect prognostic differences in BC molecular subtypes, we applied the GSVA algorithm to explore the KEGG pathways that differ across molecular subgroups. Setting |LogFc|< 0.1 and P values < 0.05 were considered as thresholds for significantly enriched pathways. Then, we conducted GO and KEGG enrichment analysis for genes that differed within molecular subtypes. Finally, we explored the differences in pathway expression between the high-and low-risk groups using GO and KEGG enrichment analysis.

Clinical characteristics of the risk scores
We scored and grouped clinicopathological characteristics for risk and applied the chi-square test or Wilcox rank sum test to analyze and compare the differential status of these clinical characteristics in different subgroups, including age, sex, TNM stage, high and low grade, and future survival status. The survival rates of different clinicopathological characteristics were then compared using K-M curves and log-rank tests.

Estimation of immune infiltration and immune checkpoint genes analysis
The proportion of different immune cell types infiltrating different BC molecular subtypes was calculated by CIBER-SORT and LM23, a novel method widely used to characterize the cellular composition of complex tissues through gene expression values in solid tumors. To reflect the ratio of immune cells to stromal cells, we also used ESTIMATE's algorithm to evaluate immune and stromal fractions. The immune cells in each sample were scored and analyzed by single-sample gene set enrichment analysis (ssGSEA) based on the above reference gene set. Then, to explore the therapeutic effect of immune checkpoint gene blockade, we also calculated the expression of key immune checkpoint genes, including PDCD1, PD-L1, and CTLA4.

Comparison of gene mutations between two risk groups
We investigated the prognostic value of high and low tumor mutation load in BC based on mutation data from the TCGA dataset using K-M analysis in the R package. The "Maftools" package in the R software was used to analyze these mutation data in both risk groups, and TMB scores were calculated for each sample. We then listed the highest mutated genes in the high-and low-risk groups separately and compared the difference in survival probability between high and low TMB. Finally, the correlation between TMB and risk scores was explored. P < 0.05 was considered statistically significant.

Construction and evaluation of the nomogram
We used the best clinically characterized variables and combined multiple risk factors to construct nomograms, which were used to estimate the overall survival of BC patients at 1, 3, and 5 years, using the "rms" R package to plot nomograms. Then, calibration curves and ROC curves were also generated to assess the stability of the models.

Screening potential drugs
Based on the risk prognostic genes associated with autophagy screened by risk model construction, we used some of the genes in the model as targets to screen potential chemotherapeutic drugs using the "pRRophetic" package in R software, and finally, identified the association between some prognostic genes and chemotherapeutic drug sensitivity.

External validation of the risk scores
To conform the accuracy of the model, we used the same model and the same coefficients in the IMvigor210 dataset as in the TCGA and 2GEO datasets, and the validation steps were the same as before we had done.

RT-qPCR and immunohistochemistry (IHC)
Total RNA was extracted from bladder cancer tissues using TRIzol reagent, and RNA concentration was calculated. cDNA was then synthesized using a cDNA reverse transcription kit. cDNA was quantified using the GAPDH housekeeping gene as the target gene reference. Gene detection was performed using SYBR Green qPCR Master Mix (Takara), and the expression of target mRNA was determined on a real-time fluorescent quantitative PCR instrument. The primer sequences of all primers are shown in Table 2. The Human Proteome Atlas (HPA, https:// www. prote inatl as. org/) has collected the entire human proteome and characterized it using immunohistochemistry and immunocytochemistry. Therefore, the protein expression of these six modeled genes between BC tissues and normal bladder tissues was analyzed using IHC from HPA database (https:// www. prote inatl as. org/).

Statistical analysis
We selected differential autophagy-associated genes by LASSO regression, univariate Cox, and multivariate Cox regression analysis, and calculated the corresponding hazard ratio (HR) and 95% confidence interval (CI). The association of risk scores with immune cells, TMB, and HLA was performed using the Spearman correlation method. Differences in overall survival  by clinicopathological characteristics were evaluated with K-M curves and then compared using log-rank test. Subsequently, ROC and calibration curves were used to assess the accuracy of the model. All statistical analyses for this study were performed using R software (version 3.6.2), IBM SPSS 20.0 software (IBM, USA), and GraphPad Prism9 software (GraphPad Software Inc., La Jolla, CA, USA), respectively. P values less than 0.05 were considered statistically significant and statistical significance was described as follows: ns: not significant; *P < 0.05; **P < 0.01; and ***P < 0.001.

Identification of differentially expressed prognostic autophagy-associated genes
All steps and processes of this study are shown in the flow chart ( Fig. 1). We identified 552 autophagy-associated BC genes in total through the dataset of BC patients in TCGA and the Human Autophagy Database (HADb). Then, according to the above autophagy-associated bladder cancer genes, using univariate Cox regression analysis and differential analysis of cancer and paracancer tissues, we obtained 117 prognosis-related genes and 148 differentially expressed genes, respectively, which were merged together and obtained 44 differential autophagy-related genes ( Fig. 2A).
The forest map showed that TRIM2, SNRNP70, SRC, IRF5, and CHMP4C were protective genes with HR values less than 1, and the rest were risk genes with HR values greater than 1 (Fig. 2B). The heat map showed the expression status of 44 genes in normal and cancerous tissues (Fig. 2C). These 44 genes were associated with each other as shown in the protein-protein interaction (PPI) network (Fig. 2D).

Construction and validation of an autophagy-related prognostic risk model
Depending on the TCGA and GSE31684/GSE15307 merged datasets with 664 samples, we finally get the expression level of 39 differential autophagy-related genes and distinguished the molecular subgroups of BC samples. Next, we made a consensus clustering analysis, and an optimal division was reached when k = 2 (C1, C2) according to the CDF curve of consensus score, with average silhouette width = 0.96. According to the curve of survival probability, C2 was better A The optimal number of clusters is calculated and divided into two clusters and their correlation area and visualization of gene clusters. B Disease-free survival, disease-specific survival, and pro-gression-free survival in C1 and C2 molecular subtypes. C Diseasefree survival, disease-specific survival, and progression-free survival in A and B gene clusters. D Distribution of cross-validation error rates and LASSO coefficients of OS-related genes than that of C1. Furthermore, the visualization of the merged gene clusters was shown (Fig. 3A). Similarly, we apply the CancerSubtypes package to classify the genes of merged datasets, and the results are shown (Supplementary document 1). We further explored the prognosis of molecular and genotypes. Molecular subtypes C1 and C2 of BC showed significant differences in disease-specific survival (DSS), progression-free survival (PFS) (log-rank P < 0.05), and no statistically significant disease-free survival (DFS) (log-rank P > 0.05) (Fig. 3B). The results for gene clusters A and B of BC, with DSS, PFS, and DFS, were similar to those of the molecular subtypes (Fig. 3C). For selecting genes associated with prognosis, we used least absolute shrinkage and choice of operator (LASSO) regression and multifactorial Cox regression algorithm (Supplementary document 2) based on the three fused datasets to construct a model (Fig. 3D). Next, the risk scores were calculated for each BC sample and the three datasets were divided into low-and highrisk groups on the basis of the median value of the risk scores. The survival status of BC patients is shown in the plot (Fig. 4A). We then compared the overall survival of patients in the whole group, training group, and test group, respectively, and found that the overall survival of patients in the low-risk group was better than that in the high-risk group (Fig. 4B). We also analyzed the models on all the datasets, the training dataset, and the test datasets separately for comparison, with the calibration curves also close to 45° (Fig. 4C). Then we conducted ROC analysis to demonstrate the sensitivity and specificity of the survival prediction model for BC patients by assessing the area under the curve (AUC). We calculated total AUCs of 0.703, 0.679, and 0.661 for 1, 3, and 5 years of follow-up, respectively; AUCs of 0.616, 0.641, and 0.619 for 1, 3, and 5 years of follow-up in the test cohort; and AUCs of 0.740, 0.697, and 0.680 in the training cohort (Fig. 4D).

Nomogram construction and validation
We explored genes with underlying differences in molecular clusters and screened 897 genes with differences between molecular subgroups using Fdr≦0.05, |LogFc|> 1 as criteria, and obtained genotyping A, B using consistency clustering analysis. We demonstrated that these 6 risk genes can equally distinguish between molecular and the genotype of BC using the 6 risk genes screened, which is consistent with the previous findings. Among the molecular subtypes, SPOCD1 and NAT8B genes were highly expressed in C2, and FKBP10, LDLR, STMN3, and ANXA2 were highly expressed in C1, and the results of gene clusters were the same as those of molecular subtypes (Fig. 5A). To establish the relationship between molecular subtypes, genotyping, risk scoring, and prognosis of BC, we risk scored each subclasses separately.
The results showed that among the molecular subtypes, the C1 typing had a higher risk score (Fig. 5B). To evaluate the relationship between patients' clinicopathological characteristics, risk scores, and overall survival, we aimed to explore whether risk scores could be used as independent predictors. The results of the univariate and multivariate Cox regression analysis (Fig. 5C) show that the risk score model we have constructed can serve as an independent prognostic factor (P value < 0.001). To predict the clinical prognosis of patients, we constructed a nomogram including age, sex, risk score, and other relevant variables to predict the survival of BC patients at 1, 3, and 5 years (Fig. 5D). The line segment of the calibration curve as shown is close to the 45° line, which suggests that the model has good predictive performance (Fig. 5E).

Functional enrichment analysis and correlations between risk score and TMB
As previously described, the survival probability of C2 subgroup of bladder cancer is superior than C1. Afterward, we obtained 83 pathways totally that lead to differences between molecules with different prognostic outcomes (Supplementary document 3), 15 pathways each of upregulated and downregulated are shown here (Fig. 6A). We observed leishmania infection, hematopoietic cell lineage, cell adhesion molecules cams, and other inflammation-related pathways were positively correlated with C1 subclass, but metabolic-related pathways such as glycerophospholipid metabolism, alpha linolenic acid metabolism, linoleic acid metabolism, and fatty acid metabolism were positively associated with C2 subgroup.
We also investigated the expression analysis of differential genes intrinsic to the molecular typing of BC and obtained 897 differentially expressed genes, and then we performed a functional enrichment analysis (Fig. 6B). In biological processes, these differentially expressed genes were mainly enriched in extracellular matrix organization, extracellular structure organization, negative regulation of immune system process, leukocyte cell adhesion, and cell chemotaxis processes. In cellular components, differential genes are mainly enriched in collagen containing, focal adhesion, cell-substrate junction. In molecular functions, they are mainly enriched in extracellular matrix structural constituent, glycosaminoglycan binding, and heparin binding, and in KEGG, cytokine-cytokine receptor interaction, P13-Akt signaling pathway, focal adhesion, proteoglycans in cancer, cell adhesion molecules, and other pathways.
Next, we further verified the richness and difference of pathway enrichment in the high-and low-risk groups and conducted GO enrichment and KEGG pathway analysis for different risk groups (Fig. 6C). GO analysis showed that extracellular matrix structural constituent, external Page 9 of 21 46 encapsulating structure organization, collagen containing extracellular matrix myeloid leukocyte migration, and granulocyte migration were enriched in the high-risk group. KEGG pathways enrichment analysis revealed that ECM receptor interactions, focal adhesion, leishmania infection, and cytokine receptor interactions were particularly significant in the high-risk group.
In addition, we explored whether TMB was associated with risk score signatures and we analyzed the relationship between risk groups and gene mutations. The status of the top 20 mutational genes was as shown (Fig. 6D). There is no difference in the type of these mutated genes, but TP53, RB1, and FGFR3 were highly variable across risk subgroups. We analyzed the relationship between risk score and  TMB using K-M survival curves, whose results showed that the prognosis was best for low risk and high TMB and worst for high risk and low TMB.

The relationship between risk scores and clinical characteristics
According to the classification of clinical characteristics, we scored the risk separately for each stratum of clinical characteristics. Based on information from the TCGA and GSE31684/GSE15307 databases, patients aged > 65 years had significantly higher risk scores in the former compared with those aged ≤ 65 years; males had higher risk scores than female patients; and high grade had higher risk scores than low-grade patients; in addition, the T-stage (T1/T2 and T3/T4) and N-stage (N0 and N1/N2/ N3) of the tumor, the former had higher risk scores than the latter (P < 0.05); regarding the future survival status The top most 10 enriched GO enrichment pathways in biological processes, cellular components, and molecular functions, respectively. C Bubble chart of GO and KEGG analysis of top 10 pathways in different risk groups. D The waterfall plot of gene mutation features established with low-and high-risk scores and K-M analysis of the OS between the low-and high-TMB groups and survival analysis among four different groups stratified by both TMB and risk score of patients, those with lower risk scores had better longterm survival status. Meanwhile, based on the molecular subclasses of BC, we performed another chi-square test for the above clinical characteristics and obtained similar results (Fig. 7A-B). We further explored the survival outcomes of patients with different clinical characteristics in different risk groups using K-M survival analysis. In these patients, a higher risk score had a lower long-term survival probability, regardless of age > 65 or ≦65 years, male or female, high grade or low grade, T-stage (T1/T2 or T3/T4), and N-stage (N0 or N1/N2/N3) of the tumor (Fig. 8). Fig. 7 A A strip chart and the scatter diagram B demonstrated that gender, tumor grade, T stage, T stages (T1/T2 versus T3/T4), and N stages (N0 versus N1/N2/N3) and future state was apparently associated with the risk score. 0 represents negative, 1-4 represent positive

The immune cell infiltration analysis, immune checkpoint genes, and immunotherapy response prediction
We used the CIBERSORT and ESTIMATE algorithms to evaluate the degree of infiltration of 23 immune cell types in different BC molecular subtypes and different risk groups. According to the results, the proportion of activated B cells, activated CD4 T cells, activated CD8 T cells, myeloid-derived suppressor cellna (MDSCna), natural killer cells (NK), neutrophils, monocytes, and regulatory T lymphocytes was significantly higher in the C1 group than in the C2 group (Fig. 9A, C). In the highrisk score group, activated CD4 T cells, activated CD8 T cells, MDSC, macrophages, mast cells, NK cells, neutrophils, plasma cells, regulatory T lymphocytes, and type.27 helper cellna were more abundant in the high-risk group than in the low-risk group. To establish the relationship between molecular subtypes, genotyping, risk scoring, and prognosis of BC, we risk scored each subclasses separately. The results showed that among the molecular subtypes, the C1 typing had a higher risk score (Fig. 9B). As shown in Fig. 9D, high-risk scores were positively associated with higher immune scores, and stromal scores. In the high-risk group, CD27, LAG3, PDCD1, CD274, CTLA4, CD276, CD44, CD48, NRP1, and PDCD1LG2 checkpoint genes were highly expressed (Fig. 9E). HLA can be classified into class I, class II, and class III gene regions, which are mainly involved in antigen presentation on the cell surface and in the immune regulation of the body. HLA can activate the immune response by delivering processed antigens to T lymphocytes. It has been shown that low expression of some molecules in HLA is associated with poor prognosis in some human cancers. As shown in the figure (Fig. 9F), the expression of HLA-A, HLA-B, HLA-C, HLA-DRA, HLA-DRB, and other genes was significantly higher in the high-risk group than in the lowrisk group. We established a correlation between immune checkpoint genes and immunophenoscore (IPS) based on the expression status of CTLA-4 and PD-1, and we compared the IPS between the high-and low-risk groups, whose results showed that PD-L1 and CTLA-4 gene positivity had significantly different immune responses in the high-and low-risk subgroups (Fig. 9G).

The sensitivity of prognostic gene expression to chemotherapy
To explore the correlation between risk genes and the efficacy of chemotherapeutic agents for BC, we used a constructed risk score model to screen these compounds for sensitive chemotherapeutic agents by analyzing the transcriptomes of 60 cancer cell lines collected from the CellMiner database based on the correlation between inhibitory concentration (IC50) and risk gene expression. All risk genes were significantly correlated with the sensitivity of some chemotherapeutic agents (P < 0.05), and the top 20 alignments with bigger |cor| are shown here (Fig. 10). Raloxifene, etoposide, camustine, arsenic trioxide, cyclophosphamide dromostanolone, tamoxifen, and teniposide were negatively correlated with increased expression of ANXA2 gene, but positively correlated with high expression of Lrofulven. For SPOCD1 gene, rapamycin, abiraterone, zoledronate, tromethamine, bleomycin, everolimus, and cabozantinib were found to be sensitive. Difttox ontak was shown to be sensitive to high expression of NAT8B gene; however, oxaliplatin was apparently insensitive to LDLR expression.

Signature validation in IMvigor210 dataset
To prove whether the autophagy-related prognostic model constructed for the 6 risk genes is reliable, we chose to perform external validation in an additional database, IMvigor210. We calculated risk scores for 348 samples in the IMvigor210 dataset using the same formula and divided them using the same cutoff values. The results showed that in the external validation set, the overall survival rate was significantly higher in the low-risk patients compared to the high-risk group of BC patients. In addition, the different risk score distributions and the comparison of survival time and status of patients in the high-and low-risk groups are shown in Fig. 11A. The ROC curves showed that the model was accurate and had good predictive performance (Fig. 11B). The expression of ANXA2, LDLR, FKBP10, and STMN3 genes was more highly expressed in the high-risk group, while the remaining genes were more highly expressed in the low-risk group, which was consistent with the results of the TCGA-2GEO fusion dataset. For immunotherapy response, the ratio of patients in complete remission (CR)/partial remission (PR) and stable disease (SD)/progressive disease (PD) in both risk groups was shown (Fig. 11C). We found a significant difference between the two groups. Subsequently and analyzed the differences in immune cell infiltration in different risk subgroups using ssGSEA, and the results showed a higher proportion of memory B cells, CD4 T cells, resting state NK cells, activated NK cells, activated mast cells, and neutrophils in the high-risk group (Fig. 11D-E). The relationships between stromal cells, immune scoring, and risk score are shown (Fig. 11F).

Validation of the expression properties of six modeled genes using RT-qPCR and HPA
The mRNA expression levels of these six modeled genes were measured in normal human bladder tissues and bladder tumor tissues by RT-qPCR. The results showed that the expression levels of four of these genes (FKBP10, LDLR, STMN3, ANXA2) were significantly different between BC and normal bladder tissues (Fig. 12A). Moreover, we compared the protein expression of model genes between normal and tumor tissues using immunohistochemical images from the HPA database, and of those, three were significantly upregulated (LDLR, FKBP10, and ANXA2) while the other 3 genes were not significantly different (Fig. 12B). Fig. 9 Immune cells infiltration, checkpoint genes, and immunotherapeutic response. A Immune cell infiltration levels between the two molecular subgroups and correlation between estimated immune cell infiltration levels and two different risk scores. B Differences in risk score between the two molecular and two gene clusters. C Spearman analysis showed the correlation between AI score and tumor immune infiltrating cells. Red represented positive correlation and blue represented negative correlation. D The relationship between risk score and stromal score, immune score. E The expression of key immune checkpoint genes between different groups. F The expression levels of HLA genes in high and low subtypes. G Correlation key immune checkpoint genes and risk score. H IPS in different group comparisons of the clinical benefit of targeting CTLA4( +) + PD-1(-), CTLA4(-) + PD-1( +), CTLA4( +) + PD-1( +), CTLA4(-) + PD-1(-) on immunotherapy in the high-or low-risk score subgroups ◂ Discussion BC is one of the most common malignancies of the urinary tract and imposes a heavy health and economic burden on society. Traditional treatment options for BC include surgery, radiotherapy, and chemotherapy, but due to the lack of a good predictive model for prognosis, the risk of death for patients with advanced BC is still high and its treatment remains a challenge. Therefore, the ability to accurately predict the prognosis of BC patients is a prerequisite for clinicians to develop personalized treatment and review plans. Autophagy plays a crucial role in cancer because it can either suppress tumorigenesis or promote tumor proliferation (Kimmelman and White 2017). The construction of autophagy-related gene models has been reported in gastric cancer (Tong et al. 2021), colon cancer (Chen et al. 2021;Zhu et al. 2020), and other tumors. However, the impact of genetic alterations in key autophagy-regulated genes on BC is unclear. Recently, new biomarkers that can be used to screen for immunotherapeutic agents have been introduced into clinical practice, which shows good therapeutic promising results to some extent, including anti-PD-1/ PD-L1 antibodies for other oncology treatments. However, the response of most BC patients to immunotherapy is still insufficiently elucidated.
In this study, we identified C1 and C2 molecular subtypes of BC based on TCGA and GSE13507/GSE31684 datasets. We found that the clinical prognosis of patients differed between these two subgroups (C2 better than C1). Then, we analyzed the genes inherent in molecular typing by consensus clustering analysis, and found that they can be divided into cluster A and B gene sets, with the different of OS, PFS, and DFS between two groups (A is better than B). Furthermore, based on the genes in the TCGA and GSE13507/GSE31684 merged datasets, we obtained 6 autophagy-related risk genes through LASSO and multivariate Cox regression analysis. Depending on the median value of risk scores, we divided the BC samples of the dataset into the high-risk and low-risk groups. K-M survival analysis showed that OS in the high-risk group was poorer in the lower-risk group. Through the above subtypes, we found that different clinical outcomes were related to their intrinsic genetic differences. In order to further study the underlying mechanism within the differences, we subsequently performed the pathway enrichment analysis of GSEA in each group. Pathways such as focal adhesion and the extracellular matrix (ECM) receptor interaction, and the low-risk group were mainly enriched for metabolicrelated pathways. It has been reported in the literature that cytokine-cytokine receptor interaction pathways play important roles in adaptive inflammatory host defense, cell growth, differentiation, cell death, angiogenesis, and developmental and repair processes aimed at restoring homeostasis (Schreiber and Walter 2010;Wang et al. 2021). In  The distribution of immune cell fraction in the low-and high-risk groups. E Correlation between risk score and tumor immune infiltrating cells. Red indicates positive relationship and blue indicates negative relationship. F Associations risk score between stromal score and immune score their relative signaling cascade activities will lead to severe damage at the cellular level of the body, which in turn will promote tumor invasion and metastasis (Guo et al. 2021;Henke et al. 2019;Schaefer and Dikic 2021). These could be potential targets for future treatments.
In recent years, the immune microenvironment has received increasing attention in cancer research, and the immune microenvironment composed of multiple immune cells plays an important role in regulating the proliferation, invasion, and migration of cancer cells Wu and Dai 2017). According to the different prognostic outcomes of bladder cancer patients in the high-risk and low-risk groups, we further analyzed the differences in immune cell infiltration between the two groups by the ssGSVA algorithm. In the high-risk group, the levels of immune activation-related cells such as CD8 + T cells, macrophages, tumor-infiltrating lymphocyte (TIL), MDSCna, mast cells, NK cells, and B cells were higher. According to the research, the M2 phenotype of macrophages is associated with tumor progression and metastasis, leading to poor prognosis of patients (Malfitano et al. 2020;Chen et al. 2018;Chanmee et al. 2014). In addition, in cancer, MDSCna promote immune evasion, stimulate proliferation, promote invasion and metastasis, and ultimately impede immunotherapy (Tran et al. 2021), and the increase of these immune cells causes changes in the immune microenvironment of tumor cells, which consists of malignant and nonmalignant cells, in which cancer cells absorb non-malignant cells, including immune cells and stromal cells, and the degree of immune cell infiltration in the TME plays an important role in regulating the proliferation, invasion, and migration of cancer cells (Siddiqui and Glauben 2022). This also explains why the prognosis of high-risk populations is poor despite the presence of a large number of immune cell infiltrates. Therefore, systematic analysis of the TME and the diversity of tumor-infiltrating immune cell species Fig. 12 Verification of the six modeled genes expression in BC and normal bladder tissue using RT-qPCR and the HPA database allows the selection of immunotherapy-sensitive patients and thus the administration of more appropriate regimens. The safety and efficacy of PD-1 immune checkpoint inhibitors and CTLA-4 inhibitors have been reported to be evaluated in several tumor treatments with promising results (Liu et al. 2021;Dong et al. 2021;Havel et al. 2019). In this study, we also evaluated the correlation between risk models and immunotherapy response in bladder cancer, and the results showed that the clinical benefits of immunotherapy against PD-1 and CTLA4 differed significantly in the lowand high-risk groups, and also indicated that the risk model we constructed could accurately reflect the therapeutic effects of immune checkpoint blockade. This also provides new directions for the mechanisms of PD-1 and CTLA-4 immune checkpoints and new therapeutic targets for bladder cancer patients with poor prognosis.
BC is a heterogeneous tumor and its treatment strategy cannot be generalized, so the identification of more effective biomarkers is necessary to personalize the treatment of BC patients (Patel et al. 2020). In this study, we analyzed the differences in TMB in two risk groups. K-M survival curves showed that tumors with high mutations in the low-risk group had a better prognosis. Theoretically, a higher TMB increases the probability of tumor neoantigen production and therefore immune recognition and tumor cell killing (Wang and Li 2019;Sholl et al. 2020). It has also been suggested that tumors with more mutated genes tend to produce more mutant RNAs and proteins that are more easily recognized by the immune system and respond well to immunotherapy (Rizvi et al. 2015). FGFR3 mutations have been reported to be common in non-invasive BC and are associated with a good BC prognosis (van Rhijn et al. 2020;Sweis et al. 2016). These mutated genes can be potential targets for therapy because of their different levels in different risk groups, so FGFR3 and RB1. HLA genes are important immune genes in the body and are mainly involved in immune responses as endogenous and exogenous antigen presenting molecules (Mcgranahan et al. 2017). Tumor-induced immune escape can alter the expression of HLA genes, allowing tumors to evade the immune system without being killed. Several reports in the literature have mentioned that HLA expression is associated with the prognosis of colorectal, gastric, and cervical cancers (Piao et al. 2021). We analyzed HLA gene expression in the high-and low-risk groups of BC patients and found that HLA gene expression was higher in the high-risk group. It has also been mentioned in the literature that high expression of HLA can affect the prognosis of BC (Kobayashi et al. 2022) and heterozygosity for HLA-I gene are related to a better prognostic outcome of immune checkpoint blockade therapy (Chowell et al. 2018. To investigate more about the advantages of our signature, we compared some other risk models simultaneously. The results showed that our model had a good application compared to other models (Lu et al. 2022;Sun et al. 2021). The model we constructed is accurate and has good predictive performance based on ROC and calibration curves, and also showed that we can predict the level of risk of patients based on the expression of these risk genes, which also provides the possibility of precise treatment for BC patients in the future. This is the first time that we have constructed a novel prognostic risk model integrating six autophagy-related genes, and more trials are needed to further validate the model in future clinical studies, as well as a larger sample and more center-based studies in the future to prove our conclusions.

Conclusions
In summary, the results of our study suggest that patients with BC can be classified into different molecular types, gene clusters, and different risk groups according to 6 autophagy-related genes. These classifications of BC patients might possibly provide more accurate personalized treatment. However, larger and multicenter studies are necessary to confirm the advantages of our model.