N6-methyladenosine-related non-coding RNAs are potential prognostic and immunotherapeutic responsiveness biomarkers for bladder cancer

Background Bladder cancer (BC) is a commonly occurring malignant tumor of the urinary system, demonstrating high global morbidity and mortality rates. BC currently lacks widely accepted biomarkers and its predictive, preventive, and personalized medicine (PPPM) is still unsatisfactory. N6-methyladenosine (m6A) modification and non-coding RNAs (ncRNAs) have been shown to be effective prognostic and immunotherapeutic responsiveness biomarkers and contribute to PPPM for various tumors. However, their role in BC remains unclear. Methods m6A-related ncRNAs (lncRNAs and miRNAs) were identified through a comprehensive analysis of TCGA, starBase, and m6A2Target databases. Using TCGA dataset (training set), univariate and least absolute shrinkage and selection operator (LASSO) regression analyses were performed to develop an m6A-related ncRNA–based prognostic risk model. Kaplan-Meier analysis of overall survival (OS) and receiver operating characteristic (ROC) curves were used to verify the prognostic evaluation power of the risk model in the GSE154261 dataset (testing set) from Gene Expression Omnibus (GEO). A nomogram containing independent prognostic factors was developed. Differences in BC clinical characteristics, m6A regulators, m6A-related ncRNAs, gene expression patterns, and differentially expressed genes (DEGs)–associated molecular networks between the high- and low-risk groups in TCGA dataset were also analyzed. Additionally, the potential applicability of the risk model in the prediction of immunotherapeutic responsiveness was evaluated based on the “IMvigor210CoreBiologies” data set. Results We identified 183 m6A-related ncRNAs, of which 14 were related to OS. LASSO regression analysis was further used to develop a prognostic risk model that included 10 m6A-related ncRNAs (BAALC-AS1, MIR324, MIR191, MIR25, AC023509.1, AL021707.1, AC026362.1, GATA2-AS1, AC012065.2, and HCP5). The risk model showed an excellent prognostic evaluation performance in both TCGA and GSE154261 datasets, with ROC curve areas under the curve (AUC) of 0.62 and 0.83, respectively. A nomogram containing 3 independent prognostic factors (risk score, age, and clinical stage) was developed and was found to demonstrate high prognostic prediction accuracy (AUC = 0.83). Moreover, the risk model could also predict BC progression. A higher risk score indicated a higher pathological grade and clinical stage. We identified 1058 DEGs between the high- and low-risk groups in TCGA dataset; these DEGs were involved in 3 molecular network systems, i.e., cellular immune response, cell adhesion, and cellular biological metabolism. Furthermore, the expression levels of 8 m6A regulators and 12 m6A-related ncRNAs were significantly different between the two groups. Finally, this risk model could be used to predict immunotherapeutic responses. Conclusion Our study is the first to explore the potential application value of m6A-related ncRNAs in BC. The m6A-related ncRNA–based risk model demonstrated excellent performance in predicting prognosis and immunotherapeutic responsiveness. Based on this model, in addition to identifying high-risk patients early to provide them with focused attention and targeted prevention, we can also select beneficiaries of immunotherapy to deliver personalized medical services. Furthermore, the m6A-related ncRNAs could elucidate the molecular mechanisms of BC and lead to a new direction for the improvement of PPPM for BC. Supplementary Information The online version contains supplementary material available at 10.1007/s13167-021-00259-w.


Introduction
Bladder cancer (BC) is a malignant tumor of the urinary system and the tenth most commonly reported cancer worldwide, demonstrating high morbidity and mortality rates [1]. BC is classified into two types depending on the existence of tumor invasion into the muscle layer of the bladder, that is, non-muscle-invasive and muscleinvasive BC types are reported. Non-muscle-invasive BC (NMIBC) accounts for approximately 70% of all newly diagnosed BC cases [2]. Transurethral resection of bladder cancer (TURBT) is the first choice of treatment, followed by intravesical bacille Calmette-Guerin (BCG) installations or chemotherapy [3]. Although the 5-year survival rate of NMIBC is approximately 90%, the postoperative recurrence rate is relatively high (50-70%) [4]. Muscleinvasive BC (MIBC) represents approximately 20% of all newly diagnosed BC cases, and exhibits a 5-year survival rate of only 50% after radical cystectomy (RC) and pelvic lymph node dissection, with or without chemotherapy and radiation therapy [2]. Although current comprehensive treatment programs such as surgery, radiotherapy, chemotherapy, and targeted therapy can help prolong the overall survival of patients to a certain extent, the overall recurrence and mortality rates of BC remain high, and the prognosis is poor [3]. As BC is highly heterogeneous, predictive, preventive, personalized medicine (PPPM) is an effective strategy used to improve treatment outcomes and patient prognosis [5]. PPPM requires the use of various effective molecular biomarkers, including early diagnostic and prognostic biomarkers that can help clinicians to identify patients in need of early, aggressive management, and predictive biomarkers that can forecast and stratify responses to emerging targeted therapies [5]. In recent years, with the application of multi-omic approaches in cancer research, many BC biomarkers have been reported [6][7][8]; however, none of them have been introduced into clinical practice. Therefore, the outcomes of PPPM for BC remain unsatisfactory.
RNA acts as a carrier as it transfers genetic information from DNA to proteins and participates in the regulation of various biological processes [9]. Similar to DNA and proteins, RNA undergoes multiple chemical modifications. Presently, more than 170 RNA modifications have been identified in all living organisms, of which methylation is the most important modification reported [10,11]. N 6 -methyladenosine (m 6 A) is widely distributed in various types of RNA, including mRNA, tRNA, rRNA, miRNA, and long non-coding RNA (lncRNA) [12,9]. It is the most commonly documented type of base modification and it is conserved in various eukaryotic organisms [11,13]. During the process of m 6 A modification, the methyl group on S-adenosylmethionine (SAM) is transferred to the sixth nitrogen atom of adenine; this process is dynamic and reversible and involves three types of regulators, demethylases (erasers), methyltransferases (writers), and methylated binding proteins (readers) [14,15]. Coupled with the action of m 6 A regulators, m 6 A modifications constitute a rapid mechanism for the coordination of RNA processing and metabolism, and for the regulation of almost all vital normal bioprocesses, including cell differentiation, embryonic development, biological rhythm regulation, heat shock response, and DNA damage repair. [16]. As expected, abnormal m 6 A modifications are closely related to the pathogenesis of multiple diseases, especially tumors [17]. Many studies have shown that abnormal m 6 A modifications, which lead to an imbalance in the expression of oncogenes and tumor suppressor genes, may contribute to the initiation and progression of tumors, and affect patient sensitivity to radiotherapy and chemotherapy, as well as clinical prognosis [10,[18][19][20].
Although non-coding RNAs (ncRNAs) do not encode proteins, they are considered important regulators of gene expression and are involved in the initiation and development of several diseases [21]. Depending on its length, RNA can be classified into two categories, namely small RNA and long-chain RNA. The former is usually less than 50 nucleotides in length, and includes tRNA, rRNA, and miRNA, while the latter is usually more than 200 nucleotides in length [22]. ncRNAs can regulate gene expression in a variety of ways and play pivotal roles in cancer development and progression [22,23]. For example, dysregulated ncRNAs possess the potential to initiate tumorigenesis, to promote invasion and metastasis, and to confer drug resistance in BC [24]. Of note, studies have demonstrated that m 6 A modification can cause ncRNA function abnormalities by regulating ncRNA processing, thereby promoting tumorigenesis [25]. For example, the methylation-binding protein, YTHDF3, can recognize the m 6 A site on lncRNA GAS5 and promote its degradation, thus activating the Hippo-YAP signaling pathway which promotes colon cancer promotion [26]. The methyltransferase, METL3, promotes the proliferation of BC cells by accelerating pri-miR221/222 maturation in an m 6 A modification-dependent manner [27]. Currently, only a few studies are available which have explored the role m 6 A ncRNA modifications in the progression of BC, and the molecular mechanisms underlying m 6 A modification in BC have not been comprehensively clarified. Therefore, understanding the role of m 6 A ncRNA modification in BC will help clarify the complex molecular mechanisms underlying BC and aid the identification of biomarkers for early diagnosis and prognostic evaluation.
In the present study, we had explored the potential application value of m 6 A-related ncRNAs in BC. Our results indicate that m 6 A-related ncRNAs may be considered novel biomarkers for predicting BC outcome and immunotherapeutic responsiveness, facilitating patients by providing them with targeted prevention and personalized medical service, thereby contributing to the improvement of PPPM for BC.

Construction and validation of the m 6 A-related ncRNA-based risk score model
Based on the TCGA dataset, through univariate Cox regression analysis, we identified 14 m 6 A-related prognostic ncRNAs. The R package "glmnet" [31] was used to perform least absolute shrinkage and selection operator (LASSO) Cox regression analysis (with the penalty parameter estimated by 10-fold cross-validation), and finally, we established a risk model for BC with 10 m 6 A-related ncRNAs. Risk score was calculated using the following formula: where n represents the number of m 6 A-related ncRNAs screened by LASSO, and a and e represent the coefficients correlated with survival and m 6 A-related ncRNA expression, respectively, as determined via LASSO.
To further evaluate the prognostic performance of this risk model, the GSE154261 dataset obtained from the GEO database was used as the testing set. Using the risk score formula, risk scores were calculated for each patient in the GSE154261 dataset, and patients were divided into high-and low-risk groups based on the median value of the prognostic risk score. The R packages "survMiner" and "survival" were used to perform Kaplan-Meier survival analysis to identify differences in OS between the two groups. The R package "pROC" was used to generate Risk score = e ∑ n i=1 a i e i the receiver operating characteristic (ROC) curve and to calculate the area under the curve (AUC) [32].

Construction and verification of the predictive nomogram
Univariate and multivariate Cox regression analyses were performed to identify independent variables, such as age at diagnosis, sex, race, pathological grade, risk score, and clinical stage. Furthermore, to individualize the predicted 1-and 2-year survival probabilities, the R package "rms" was used to generate a nomogram that included significant clinical characteristics and calibration plots. Correction curves based on the Hosmer-Lemeshow test were generated to compare prediction accuracy between the observed and model-predicted outcomes.

Functional enrichment analysis of differentially expressed genes (DEGs)
DEGs between the high-and low-risk groups were identified based on the conditions, | log 2 (fold change)| > 1 and FDR < 0.05, using the R package "DESeq2" [33]. DAVID (version 6.8, https:// david. ncifc rf. gov/) was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment and cluster analyses [34].

Evaluation of the immunotherapeutic responsiveness predictive ability of the risk model
We calculated the risk scores of 348 patients in the "IMvig-or210CoreBiologies" dataset. Based on the median values, patients were divided into high-and low-risk groups. Survival analysis was performed to identify differences in the survival status between the two groups. Additionally, we assessed the relationship between risk score and immunotherapeutic responsiveness.

Identification of m 6 A-related ncRNAs in BC patients
The study flowchart is illustrated in Fig. 1a. First, we identified m 6 A-related ncRNAs by comprehensively analyzing the results of two databases (starBase and m 6 A2Target) and by performing correlation analyses on m 6 A regulators and ncRNAs (miRNAs and lncRNAs) using the TCGA dataset. After performing univariate Cox regression analysis and LASSO regression analysis, an m 6 A-related ncRNA-based risk score model was constructed. We verified the prognostic evaluation performance of this risk score model using TCGA (training set) and GEO (testing set) datasets. Finally, based on the risk model established herein, patients in the TCGA dataset were divided into the high-and low-risk groups. We further analyzed differences in the survival status, m 6 A regulators, gene expression patterns, and immunotherapeutic responsiveness between the two groups. Through Pearson's correlation analysis of ncRNAs and m 6 A regulators in the TCGA dataset, we obtained 13,596 lncRNAs and 1306 miRNAs that satisfied the filtering criteria. Figure 1b and 1c depict the top 20 miRNAs and lncR-NAs, respectively, that were most significantly associated with m 6 A regulators in the TCGA dataset. Overall, lncRNAs were more significantly correlated with m 6 A regulators than miRNAs. Among them, MIR647 and lncRNA OTUD6B-AS1 were found to be significantly correlated with m 6 A regulators METTL3 and YTHDF3, with Pearson's correlation coefficients of 0.60 (P < 0.001) and 0.61 (P < 0.001), respectively. Detailed results of the correlation analysis between ncRNAs and m 6 A regulators are listed in Supplementary Tables 1 and 2. Using the m 6 A2Target database, a total of 704 potential downstream m 6 A regulator target genes were identified, including 417 ncRNAs. Using the starBase database, we identified 3778 ncRNAs that established interactions with m 6 A regulators. Then, ncRNAs obtained from three datasets (i.e., the m 6 A regulator and ncRNA correlation analysis, starBase, and m 6 A2Target databases) were investigated and they were found to overlap, and 183 unique m 6 A-related ncRNAs were identified (Fig. 2a).
Next, we assessed the prognostic evaluation power of the risk model in the TCGA dataset. Based on the median value of the risk score (0.934), samples were divided into high-risk (n = 203) and low-risk (n = 204) groups. Figure 3a depicts the risk score distribution between the patients. Survival analysis showed that there were significant differences in the survival status between the two groups (P < 0.01), and the survival time of the low-risk group was longer than that of the high-risk group (Fig. 3b). Assessment of the predictive power of the risk score on survival showed that the maximum AUC of the ROC curve was 0.62 (Fig. 3c). Additionally, ROC analysis revealed that the risk model exhibited a good predictive power for the probability of OS at 3 and 5 years for the TCGA dataset (Fig. 3c).

Validation of the prognostic evaluation power of the m 6 A-related ncRNA-based risk model using the testing set
To further explore the prognostic evaluation efficiency of this risk model, we adopted the same method as described above to analyze the GSE154261 dataset (testing set) obtained from the GEO database. Figure 3d depicts the risk score distribution of 73 samples in the GSE154261 dataset. Based on the median value of the risk score (0.514), the samples were divided into high-risk (n = 36) and low-risk (n = 37) groups. As shown in Fig. 3e, there was a significant to m 6 A regulators in TCGA. c Heatmap for the top 20 lncRNAs that were most related to m 6 A regulators in TCGA. DEGs: Differentially expressed genes difference in survival time between the two groups, with the low-risk group exhibiting a longer survival time than the high-risk group (P = 0.03). This finding was consistent with that of the survival analysis performed using the TCGA dataset. Assessment of the predictive power of the risk score on survival showed that the maximum AUC of the ROC curve was 0.83 (Fig. 3f). The risk model showed excellent predictive power for the probability of OS at 1, 3, and 5 years (Fig. 3f).

The m 6 A-related ncRNA-based risk model is an independent prognostic factor for BC
To determine whether the risk model was an independent prognostic factor, the risk score was considered as a new variable and we performed univariate and multivariate Cox regression analyses; other clinical characteristics such as age, sex, pathological grade, race, and clinical stage were also considered. The results of the univariate Cox regression analysis showed that the risk score (P < 0.001, hazard ratio [HR] = 0.58, 95% CI: 0.43-0.78), as well as age (P = 0.001, HR = 0.52, 95% confidence interval [CI]: 0.35-0.77) and clinical stage (P < 0.001, HR = 0.46, 95% CI: 0.32-0.66), was significantly correlated with OS (Fig. 4a). Additionally, a higher risk score indicated poorer OS in BC patients. Multivariate Cox regression analysis results further confirmed that a "high risk score" was an independent factor associated with poor OS (P = 0.003, HR = 0.63, 95% CI: 0.47-0.86) (Fig. 4b).

Construction and validation of the prognostic nomogram
To provide clinicians with a quantitative approach for predicting the prognosis of BC, a nomogram that was constructed by integrating risk score and other independent clinical risk factors (age and clinical stage) was established and used for further analysis. The line segment corresponding to each independent risk factor in the nomogram is marked with a scale, which represents the range of each independent risk factor, and the length of the line segment reflects the contribution of each independent risk factor. The value of each independent risk factor corresponding to the first row of "Points" is the single score. The total score obtained by adding up single scores will eventually correspond to the "Total Points" row. Finally, according to the position of total score at different observation times, OS is predicted. Poor prognosis was represented by a higher total number of points on the nomogram (Fig. 5a). We also compared the predictive accuracy of this nomogram with that of the analysis based on age and risk score and found that the performance of the nomogram (AUC = 0.83) was better than the analysis conducted using data on age (AUC = 0.62) and risk score (AUC = 0.69) (Fig. 5b). Additionally, predictions proposed by using the calibration curve of the nomogram for the first-, second-, and third-year OS were closer to the observed OS (Fig. 5c-e).

Evaluation of the prognostic risk model and the clinical characteristics of BC
Having confirmed that the risk model was an efficient tool for predicting prognosis, we aimed to ascertain if it could reflect the clinical characteristics of BC. To achieve this purpose, we analyzed the correlation between risk score and the clinical characteristics of BC. As shown in Fig. 6a, patients older than 60 years of age presented with higher risk scores than those younger than 60 years; female patients presented with higher risk scores than male patients; Black or African American (BAA) and White patients presented with higher risk scores than Asians; patients with higher BC pathological grades presented with higher risk scores than those with lower BC pathological grades; as the tumor TNM grade increased, the risk score also tended to increase; patients who presented with no initial treatment response (NR) exhibited higher risk scores than those who demonstrated an initial treatment response. These results indicated that the risk model was related to the progression of BC. Subsequently, we conducted a survival analysis for the OS of subgroups with different clinical characteristics. We found that the OS of patients in the low-risk group was better than that of patients in the high-risk group, in the "older than 60 years" group, in the "therapy response" group, in the male patients' group, and in patients in three races (BAA, White, and Asian) (Fig. 6b).

Evaluation of molecular network alterations between the high-and low-risk groups
Owing to the significant differences in survival status observed between the high-and low-risk groups, we analyzed the gene expression profiles of the two groups. According to the previously established criteria for DEGs, we obtained data on 1058 DEGs (Supplementary Table 3 and Fig. 3a). As shown in supplementary Fig. 3b, DEGs exhibited significantly different expression patterns between the high-and low-risk groups. To understand the molecular network alterations involved in these DEGs, we performed cluster analysis using the results of the GO and KEGG enrichment analyses. Each cluster represented a type of molecular network system with similar functions during carcinogenesis. We found that the 3 clusters were significantly enriched ( Table 1). Cluster 1 was significantly associated with the cellular immune response, which comprised multiple biological processes and pathways related to the activation and chemotaxis of immune cells, such as neutrophil chemotaxis, monocyte chemotaxis, lymphocyte chemotaxis, and cellular response to interferon-gamma and interleukin-1. Cluster 2 was associated with cell adhesion and the interaction between cells and the extracellular matrix (ECM). Cluster 3 was involved in cellular biological metabolism, including the metabolic processes of a variety of important compounds, such as steroids, retinol, linoleic acid, and arachidonic acid. Additionally, the cytochrome P450 pathway, which is involved in drug metabolism, belonged to this cluster. m 6 A regulator and m 6 A-related ncRNA expression differences between the high-and low-risk groups m 6 A regulators, which are critical molecules in the modification process, also showed expression differences between the high-and low-risk groups. We discovered 8 m 6 A regulators; of these, IGF2BP2, IGF2BP3, and ALKBH5 were highly expressed in the high-risk group, while YTHDF1, METTL3, YTHDF2, YTHDC1, and ZCCHC4 were highly expressed in the low-risk group (Fig. 7). Additionally, among the 14 OS-related ncRNAs, 12 (MIR324, MIR93, MIR331, MIR25, AL022311.1, AC023509.1, AC012615.1, AL021707.1, AC026362.1, GATA2−AS1, AC012065.2, and HCP5) were differentially expressed between the high-and lowrisk groups. Apart from HCP5, which was highly expressed in the high-risk group, the other ncRNAs were not highly expressed in this group (Supplementary Fig. 4).

The risk model may aid the selection of patients who may benefit from immunotherapy
Immunotherapy has revolutionized the treatment of BC; however, not all patients can derive benefits. Screening for potential beneficiaries of BC immunotherapy is an effective means of improving disease prognosis and achieving precise medical treatment for BC. Hence, we evaluated the ability of the risk model to be used to predict BC immunotherapeutic responsiveness. As described previously, the risk score formula was used to calculate the risk score for each sample (n = 348) in the BC immunotherapy dataset, and based on the median value of the risk score (0.94), the samples were divided into the high-risk (n = 174) and lowrisk (n = 174) groups. Survival analysis revealed that there were significant differences (P = 0.03) in survival status between the two groups, with the survival time of patients in the high-risk group being shorter than that of patients in the low-risk group (Fig. 8a). Concerning immunotherapeutic effects, the high-risk group presented with a lower response rate (RR, the ratio of patient with complete response (CR) and partial response (PR)) than the low-risk group (20.81% vs. 24.83%) (Fig. 8b). The above result indicated that there were differences in immunotherapy responsiveness between the high-and low-risk groups. Furthermore, patients with disease progression (PD) presented with higher risk scores than those observed in CR, PR, or those with stable disease (SD) (Fig. 8c). In summary, our results indicated that this risk model could be a promising biomarker to predict immunotherapeutic responses.

Discussion
The advent of multi-omics technology has advanced our understanding of the molecular mechanisms of tumorigenesis in a profound manner, and has resulted in changes in the clinical treatment strategies adopted for several tumors [6]. Few molecular biomarkers have been established as an indispensable approach of PPPM. For example, PSA and AFP have been used as conventional biomarkers for the diagnosis of prostate and liver cancer, respectively [35,36]. EGFR serves as a biomarker for tumor pathological classification and targeted therapy in lung adenocarcinoma [37]. Research conducted in the last 30 years, which has been focused on the molecular mechanisms of BC, has provided us with information on a substantial number of potential molecular biomarkers, including early diagnostic biomarkers (such as TERT promoter mutations and chromatin-modifying gene alterations) [7], prognostic biomarkers that can be used to identify high-risk patients who require active surveillance and early aggressive treatment (such as p53, pRB, p21, p27, cyclins D1 and D3, fibroblast growth factor receptor 3 (FGFR3), and Ki-67) [38], and predictive biomarkers which are used to forecast and stratify the patient's response to chemotherapy or targeted therapy (such as EGFR, VEGFR, Bcl-2, EMMPRIN, and survivin) [38]. In addition to the Fig. 7 Expression pattern differences of m 6 A regulators between the high-and low-risk groups in TCGA above-mentioned BC tissue-based molecular biomarkers, a wide range of potential molecular biomarkers have been discovered and reported in urine [39], circulating tumor cells [40], and exosomes [41]. Unfortunately, due to the lack of availability of sufficient validation and prospective studies, no molecular biomarker has been ratified for use in the clinical management of BC. Presently, pathological grade remains a decisive indicator for evaluating the prognosis of BC. However, in clinical practice, it is not uncommon for patients with the same pathological grade to exhibit significant differences in prognosis. This implies that due to the heterogeneity of BC, the current prognosis assessment system cannot meet the needs of personalized medicine. Therefore, for BC, which demonstrates high recurrence and mortality rates with limited diagnostic and treatment methods, the identification of novel biomarkers to improve PPPM is the need of the hour. Previous studies have shown that both m 6 A modification and ncRNAs are involved in the initiation and development of a variety of tumors [42], and can therefore be considered as potential molecular biomarkers [16,22]. Of note, there exists a complex interactive regulation between m 6 A modification and ncRNAs. Studies have shown that m 6 A modification can promote tumor cell proliferation, invasion, and metastasis by regulating the splicing, maturation, transport, and stability of ncRNAs. For example, as a "writer," METTL3 can not only promote the progression of pancreatic cancer by regulating the maturation of miR-25-3P [43], but it can also promote the tumorigenesis and metastasis of nasopharyngeal carcinoma by enhancing the stability of lncRNA FAM255A [44]. Furthermore, ncRNAs may target m 6 A regulators as competitive endogenous RNAs, thereby affecting tumor progression.
For example, low miR-1266 expression promotes colorectal cancer progression by regulating demethylase FTO [45]. LncRNA GAS5-AS1 establishes interaction with demethylase ALKBH5 to increase the stability of GAS5, thereby suppressing the growth and metastasis of cervical cancer [46]. Therefore, we believe that focusing research on m 6 A-related ncRNAs may provide new opportunities for the identification of various molecular biomarkers and for the development of targeted drugs. Recent studies have shown that m 6 A-related lncRNAs can be used as prognostic biomarkers in lower-grade gliomas and lung adenocarcinomas [47]. However, studies based on the exploration of the molecular mechanisms and clinical applications of m 6 A-related ncRNAs in BC are few. Therefore, we attempted to use bioinformatics methods to construct the first m 6 A-related ncRNA-based prognostic risk model for BC.
In this study, we identified 183 m 6 A-related ncRNAs, of which 14 were associated with OS in BC. Finally, a prognostic risk assessment model consisting of 10 m 6 A-related ncRNAs was constructed. This model showed good prognostic evaluation performance in both the training and testing sets, and multivariate Cox regression analysis results revealed that the risk model could be deemed an independent prognostic factor for OS. The nomogram integrated independent risk factors that showed adequate prognostic evaluation performance and perfect consistency between the observed and predicted rates for the 1-year-, 3-year-, and 5-year OS. The above results demonstrated that the m 6 A-related ncRNA-based risk model had excellent performance in predicting patient outcomes and can be used as a potential prognostic biomarker. With its aid, we can identify high-risk patients early and implement targeted prevention through Differences in immunotherapy response between the high-and low-risk groups. c Distribution of risk score in different immunotherapy response groups. CR: complete response, PR: partial response, SD: stable disease, PD: progressive disease more active surveillance and aggressive management. For example, in NMIBC patients, we can increase the frequency of cystoscopy and extend the duration of bladder intravesical instillation. In addition, for MIBC patients, we can expand the scope of lymph node dissection and combine chemotherapy and radiotherapy after surgery. The advent of immune checkpoint inhibitors (ICIs) has revolutionized the treatment of BC, thereby ushering hope of a cure to patients with advanced BC [2]. However, only 20-30% of the patients with advanced BC will have long-lasting clinical benefits when subjected to treatment with ICIs [48]. Screening for potential beneficiaries of BC immunotherapy is an effective means of improving disease prognosis and achieving personalized medicine for BC. Of note, the risk model could be used to predict the outcomes of immunotherapy. Patients in the low-risk group were more likely to benefit from immunotherapy. Based on this model, we can preliminarily screen out patients who may benefit from immunotherapy, thereby reducing the economic burden on patients and society. Simultaneously, personalized medicine for BC patients can be achieved. Therefore, it is apparent that m 6 A-related ncRNAs are promising molecular biomarkers and can contribute to PPPM for BC.
Most m 6 A-related ncRNAs in risk models are closely related to cancer progression; for example, lncRNA HCP5 promotes human BC cell invasion and migration by targeting miR-29b-3p [49]; lncRNA GATA2-AS1 has been identified as a colon adenocarcinoma-related lncRNA [50] and is also involved in the growth regulation of non-small cell lung cancer [51]; lncRNA BAALC-AS1 regulates the proliferation of esophageal squamous cells by participating in the lncRNA BAALC-AS1/G3BP2/c-Myc feedback loop [52]; another study indicates that miR191 is abnormally expressed in more than 20 different malignancies [53], and that overexpressed miR191 in NMIBC can cause a significant decrease in EGR1 levels [54]. However, apart from lncRNA HCP5, much remains unknown regarding the underlying molecular mechanisms of the other 9 m 6 A-related ncRNAs involved in the progression of BC, and further studies are warranted.
We observed that the 1058 DEGs between the high-and low-risk groups were mainly enriched in molecular network systems associated with cellular immune responses (such as the chemokine-mediated signaling pathway, cellular response to interleukin-1, neutrophil chemotaxis, monocyte chemotaxis, and lymphocyte chemotaxis), and cellular metabolism (e.g., drug and xenobiotic metabolism by cytochrome P450, and steroid metabolism). On the one hand, the enrichment results explained and highlighted, to a certain extent, the potential applicability of the risk model in the prediction of immunotherapeutic responsiveness. On the other hand, they indicated that the biological processes related to the cellular immune response and biological metabolism should be the focus of future research to identify the role played by ncRNAs and their interactions established with m 6 A-related genes in BC.
This study presents with certain limitations. For example, although we used a testing set to evaluate the risk score model, there were fewer samples in this testing set, and our study was a retrospective study. In order to solve this limitation, we have reached agreements with multiple urology research centers in Guangdong Province to collect BC patient clinical information and sequencing data to build our own database, which can be used as a validation data set to evaluate research results based on public databases. Additionally, the specific role of m 6 A-related ncRNAs in BC risk models remains unclear, and further research should be conducted. We hope that through the publication of this paper, more researchers can pay attention to m 6 A-related ncRNAs and contribute to revealing the role of m 6 A-related ncRNAs in the occurrence and development of BC.

Conclusions and expert recommendations
In conclusion, our study is the first to explore the potential application value of m 6 A-related ncRNAs in BC. We showed that the m 6 A-related ncRNA-based risk model demonstrated excellent performance in determining patient prognosis and immunotherapeutic responsiveness. Additionally, this model can aid in providing patients with targeted prevention and personalized medical services, thus projecting a new direction to improve PPPM for BC.

Expert recommendations and outlook
The importance of sequencing technology for life science research is self-evident. At present, genomic and transcriptomic analysis of tumor tissues obtained by intraoperative resection and puncture biopsy has become a routine diagnosis and treatment of a variety of tumors. Previous studies based on sequencing data had paid more attention to various ncRNAs and m 6 A regulators, which had been proved to play an important regulatory role in carcinogenesis. This study is the first to explore the potential application value of m 6 A-related ncRNAs in BC, and we recommend this article to emphasize the importance of m 6 A-related ncR-NAs in the basic and translational research for PPPM in BC.
Sequencing data for m 6 A-related ncRNAs analysis is not limited to tumor tissues sources. For urinary system tumors, especially BC, we think that exfoliated tumor cells (ETCs) in urine could provide a new direction for the study of m 6 A-related ncRNAs in BC. The research based on ETCs has the advantages of absolute non-invasive and high patient compliance, and shows a broad clinical application prospect.
Availability of data and material Data and download URLs involved in this study had been described in detail in the materials and methods section.
Code availability Bioinformatics analysis of this study was based on R 3.5.1, and the involved R packages were described in detail in the materials and methods section. All codes were written and reviewed by M.L. and X.Z., respectively. M.L. and X.Z. had the ultimate power of interpretation and ownership of the code.
Author contribution M.L. was responsible for bioinformatics analysis, prepared figures and tables, and designed and wrote the manuscript. H.Z., B.L., D.L., and W.L. carried out data download and preprocessing. D.L. also provided some basic code and revised this manuscript. X.C. proofread the manuscript and figures. X.Z. conceived the concept, instructed bioinformatics analysis, supervised results, and was responsible for its financial supports and the corresponding works. All authors approved the final manuscript.

Declarations
Ethical approval Since the sequenced data generated from TCGA and GEO were all publicly available data, further approval by an ethics committee was not required.
Consent to participate All authors voluntarily participated in this study.

Consent for publication
All authors agreed to the publication of this paper.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.