Integrative machine learning algorithms for developing a consensus RNA modi ﬁ cation-based signature for guiding clinical decision-making in bladder cancer

Objectives: Dysregulation of RNA modi ﬁ cations has emerged as a contributor to cancer, but the clinical implication of RNA modi ﬁ cation-related genes remains largely unclear. The study focused on well-studied RNA modi ﬁ cation modalities (m 6 A, m 1 A, m 5 C and m 7 G) in bladder cancer, and proposed a machine learning-based integrative approach for establishing a consensus RNA modi ﬁ cation-based signature. Methods: Multiple publicly available bladder cancer cohorts were enrolled. A novel RNA modi ﬁ cation-based classi ﬁ cation was proposed via consensus clustering analysis. RNA modi ﬁ cation-related genes were subsequently selected through WGCNA. A machine learning-based integrative framework was implemented for constructing a consensus RNA modi ﬁ cation-based signature. Results: Most RNA modi ﬁ ers were dysregulated in bladder tumours at the multi-omics levels. Two RNA modi ﬁ cation clusters were identi ﬁ ed, with diverse prognostic outcomes. A consensus RNA modi ﬁ cation-based signature was established, which displayed stable and powerful e ﬃ cacy in prognosis estimation. Notably, the signature was superior to conventional clinical indicators. High-risk tumours presented the activation of tumourigenic pathways, with the activation of metabolism pathways in low-risk tumours. The low-risk group was more sensitive to immune-checkpoint blockade, with the higher sensitivity of the high-risk group to cisplatin and paclitaxel. Genes in the signature: AKR1B1 , ANXA1 , CCNL2 , OAS1 , PTPN6 , SPINK1 and TNFRSF14 were specially expressed in distinct T lymphocytes of bladder tumours at the single-cell level, potentially participating in T cell-mediated antitumour immunity. They were transcriptionally and post-transcriptionally modulated, and might become potentially actionable therapeutic targets. Conclusions: Altogether, the consensus RNA modi ﬁ cation-based signature may act as a reliable and hopeful tool for improving clinical decision-making for individual bladder cancer patients.


Introduction
Bladder cancer is among the most commonly diagnosed malignant tumours of the urinary system [1].It ranks higher in men than in women and is the sixth most common cancer and the ninth leading cause of cancer-related deaths among men [2].Bladder cancer has many associated risk factors, among which cigarette smoking is the most common contributor to the increased incidence of bladder cancer in Western countries [3].Bladder cancers are divided into non-muscle-invasive and muscle-invasive types with heterogeneous clinical outcomes.Patients with non-muscleinvasive cancers are mostly curable, and the 5-year overall survival (OS) rate is nearly 90 %, whereas the 5-year OS rate for patients with muscle-invasive cancers is approximately 60-70 % [4,5].Tumour-node-metastasis staging instituted by the American Joint Committee on Cancer (AJCC) has a relatively high prognostic performance and acts as the basis for clinical decision-making; however, it is insufficient to cover the clinical characteristics of all patients [6].
Surgical resection, radiotherapy, and chemotherapy are the primary treatments for bladder cancer.However, not all patients benefit from these treatments; some may experience tumour recurrence or metastasis because of bladder cancer cell muscle invasion and drug resistance [3].Clinically, immunotherapeutics, such as adoptive T cell transfer and immune checkpoint blockade (ICB), can control cancer through antigen-specific T cells [7].Owing to off-target effects and the limitations of a single target, the clinical outcomes of immunotherapy approaches are far from satisfactory [8].Bladder cancer is a complex disease with both inter-and intra-tumour heterogeneity, the clinical outcomes of which may differ greatly among patients with the same stage and treatment strategies [9].Therefore, it is important to ensure individualised treatment by identifying reliable biomarkers to optimise the prognosis and benefit of drug therapies for bladder cancer.
The N 6 -methyladenosine (m 6 A) modification of mammalian mRNAs was discovered in the 1970s and was identified as the most abundant internal modification in mRNA [10].Subsequently, with the rapid development of high-throughput technologies, research has focused on epigenetic modifications.m 6 A influences many biological processes, not only at the molecular level, such as pre-mRNA splicing, mRNA degradation, and protein translation initiation [11,12], but also at the individual level, such as development, metabolism, sex differentiation, and tumourigenesis [13,14].With the discovery of the m 6 A demethylase fat mass-and obesity-associated protein, studies on other RNA methylation modifications, such as N 1 -methyladenosine (m 1 A), C 5 -methylcytosine (m 5 C), and N 7 -methylguanosine (m 7 G), have emerged [15,16].Recent studies [17][18][19][20] have shown that modifications of m 6 A, m 1 A, m 5 C, and m 7 G are involved in the initiation and development of bladder cancer.Hence, the dysregulation of these RNA modifications may contribute to bladder cancer development.
In recent years, deep-learning-based models have emerged and gained popularity among researchers owing to their ability to automatically extract high-order information [21].With the development of bioinformatics technology, a multitude of prognostic gene signatures, particularly those derived from messenger RNA, long non-coding RNA, and microRNAs (miRNA), have been discovered and validated as potential biomarkers for bladder cancer [22][23][24].In this study, we propose a machine learning-based integrative approach to establish a consensus RNA modification-based signature in patients with bladder cancer from independent public datasets.This study may help optimise personalised treatment and improve the clinical prognosis of patients with bladder cancer.

Interaction analysis of RNA modifiers
The relationships between RNA modifiers at the mRNA level were analysed using the Pearson correlation method.RNA modifiers were imported into the STRING website (http://string-db.org/)[29].Subsequently, a protein-protein interaction network combined with enriched pathways was visualised.

Consensus clustering analysis
Based on the expression profiles of the prognostic RNA modifiers, a resampling-based consensus clustering approach was performed on the TCGA-BLCA cohort using the ConsensusClusterPlus package (http:// www.bioconductor.org/)[33].The proportion of ambiguously clustered (PAC) pairs, tracking plot, principal component analysis (PCA), and consensus heatmap were used to estimate the optimal cluster number (k).

Weighted correlation network analysis (WGCNA)
A co-expression network of TCGA-BLCA samples was constructed using the WGCNA package (http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA) [34].The appropriate softthreshold power was determined to satisfy the criteria for a scale-free network.A weighted adjacency matrix was transformed into a topological overlap matrix, followed by the generation of the corresponding dissimilarity.A dynamic tree-cutting method was adopted to merge the co-expression modules.The module with the strongest association in the Pearson correlation analysis was chosen as the RNA modificationrelated module.Genes in the module with gene significance >0.4 and module membership >0.6 were defined as RNA modification-related genes.

Integrative machine learning algorithms for a consensus signature
To develop a consensus signature with high accuracy and stable efficacy, we integrated multiple machine-learning algorithms and their combinations, as previously described [35].Univariate Cox regression was used to select prognostic RNA modification-related genes in TCGA-BLCA samples (p<0.05).Subsequently, machine learning algorithms and their combinations were implemented to fit the prediction signatures using the leave-one-out cross-validation framework in TCGA-BLCA samples.Each signature was then analysed in the four validation cohorts, and the concordance (C)-index values for each signature as well as the mean C-index of the whole cohort were computed.The signature with the highest mean C-index was considered optimal.

Immunotherapy response estimation
Homologous recombination deficiency (HRD), which is a common molecular feature of genomic instability [46], tumour mutational burden (TMB), which refers to the number of somatic mutations per megabase of genome sequences [47], single nucleotide variant (SNV) neoantigens [48], and Tumour Immune Dysfunction and Exclusion (TIDE) [49] were adopted to estimate ICB responsiveness.The cancer-immunity cycle contains a series of steps required for the immune-based control of tumour overgrowth [50].The activity at each step was computed using the GSVA package (http://biocc.hrbmu.edu.cn/TIP/)[37].

Sensitivity analysis of chemotherapy drugs
By performing pRRophetic analysis [51], the half-maximal inhibitory concentration (IC 50 ) of the chemotherapeutic drugs was calculated based on the Genomics of Drug Sensitivity in Cancer cell line expression spectra [52].

Single-cell transcriptome analysis
Single-cell RNA sequencing profiles of T cells from human bladder tumours were acquired from GSE149652 (n=28) [53].After quality control, single cells with >20 % mitochondrial unique molecular identifiers were removed using the Seurat package (https://cran.r-project.org/web/packages/Seurat/index.html)[54].The first 1,500 genes with high variability were selected for further analysis.T cell clusters were identified using the FindClusters function, which was mapped onto UMAP.Next, T cell types were annotated based on known cell markers from the CellMarker database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/)[55].

Statistical analysis
All of the analyses were conducted using the appropriate R package (version 4.2.1, https://www.r-project.org/).Continuous data between two groups were compared using the Student's t-test or Wilcoxon rank-sum test.Pearson's or Spearman's test was used for correlation analysis.Survival analysis was conducted using the univariate Cox regression or Kaplan-Meier method.The area under the curve (AUC) from the receiver operating characteristic (ROC) curve and C-index were calculated.Statistical significance was set at p<0.05.

Multi-omics landscape of RNA methylation modifiers among bladder tumours
To elucidate the features and implications of RNA modifications in bladder cancer, we conducted a multi-omics analysis of the RNA modifiers m 6 A, m 1 A, m 5 C, and m 7 G.Most of the modifiers displayed remarkable differences in transcriptional expression between normal and bladder cancer tissues (Figure 1A).Several RNA modifiers occurred frequently in somatic mutations in bladder cancer, including TET1 (22.6 %), EIF4G1 (15.5 %), ALYREF (13.1 %), LRPPRC (10.7 %), and GEMIN5 (10.7 %) (Figure 1B).Widespread CNVs of RNA modifiers were also detected (Figure 1C and D).For instance, EIF4G1 and EIF4E2 showed the highest copy number amplification and deletion frequencies, respectively.We investigated the notable interactions between RNA modifiers at both the mRNA and protein levels, for example, the LRPPRC and GEMIN5 relationship and the CYFIP1 and EIF4G1 relationship (Figure 1E and F).Additionally, RNA modifiers participated in the p53 signalling pathway, spliceosome formation, and RNA transport, revealing their importance in bladder cancer.

RNA modification-based consensus clusters in bladder cancer
Based on the univariate Cox regression results of the transcriptomes of prognostic RNA methylation modification factors (p<0.05),we combined the SNP and CNV data of corresponding genes (Figure 1G) to perform a consistent clustering analysis of multi-omics data, and divided the bladder cancer samples into categories 2-9.Consensus clustering analysis was conducted with a similarity threshold of 97 %.Two clusters had the lowest PAC values, reflecting the near-perfect stable partitioning of the samples at k=2 (Figure 1H).Moreover, the tracking plot, PCA, and consensus heatmap displayed a relatively stable classification of the samples at k=2 (Figure 1I-K).Remarkable transcriptional heterogeneity in prognostic RNA methylation modifiers was observed in the two consensus clusters (C1 and C2), with mostly higher expression in C2 than in C1 (Figure 1L).Hence, C2 tumours have more abundant RNA modifications than C1 tumours.Additionally, the heterogeneity in prognosis was evaluated in the two clusters.As illustrated in Figure 1M, the OS outcomes were poorer in C2 than in C1.

Establishment of an RNA modificationrelated consensus signature
Before conducting the WGCNA analysis, we used a clustering method to remove outlier sample points.We first selected a cut tree height to remove the outlier sample points, then set the cut tree height to 120,000 and used the branches below this cut tree height.WGCNA was performed to build a co-expression network (Figure 2A).The soft-threshold power was selected as 15 (scale-free R 2 =0.85) to meet scale-free network standards (Figure 2B).Next, the five co-expression modules were merged (Figure 2C).Correlation analysis of the modules with RNA modification classification was performed.The turquoise module showed the strongest association with the RNA modification classification (r=0.66 and p=9e-53) (Figure 2D).In this module, the correlation coefficient of gene significance with module membership was 0.86, suggesting the excellent quality of the RNA modification-related co-expression module establishment.Genes in the turquoise module with gene significance >0.4 and module membership >0.6 were regarded as RNA modificationrelated genes (Figure 2E).Among them, 176 genes were linked to patient prognosis (Supplementary Table 3), which were then subjected to a machine learning integrative procedure to construct a consensus signature.In this study, 74 predicted signatures were fitted using a leave-one-out cross-validation framework.Next, the C-index of each signature was computed for the four validation cohorts: GSE13507, GSE32548, GSE48075, and GSE69795 (Figure 2F).To screen the important features of the random forest algorithm, we used a multi-factor Cox regression model to construct a prognostic model.Intriguingly, the best prediction signature was random survival forest (RSF) with seven feature variables (PTPN6, AKR1B1, TNFRSF14, OAS1, ANXA1, SPINK1, and CCNL2), achieving the highest C-index of 0.656 in mean values of the four cohorts.

The RNA modification-related consensus signature is superior to clinical variables in estimating prognosis
The risk score for each patient with bladder cancer was computed based on the expression of RSF-derived feature variables, weighted by the corresponding regression coefficients acquired from the Cox model.Subsequently, patients were classified into low-or high-risk groups according to the median risk score.Low-risk patients had more favourable OS outcomes than high-risk patients in the TCGA-BLCA cohort (Figure 3A), which was externally confirmed in the four validation cohorts (Figure 3B-E).ROC curves were plotted to assess the discrimination of the consensus signature, with one-, three-, and five-year AUC values exceeding 0.90 in the TCGA-BLCA cohort (Figure 3F-H).Discrimination was confirmed in the validation cohorts.The C-index was >0.9 in the TCGA-BLCA cohort, and >0.6 in all validation cohorts (Figure 3I).These indicators demonstrated that the consensus signature had stable and reliable efficacy in prognostication.Furthermore, we compared the efficacy of the consensus signature with conventional clinical parameters (age, sex, stage, and grade).The detailed clinical features are displayed in Table 1.The consensus signature presented notably superior accuracy in prognosis estimation (Figure 3J and K).

Heterogeneity in genetic alterations, and signalling pathways between low-and high-risk tumours
Bladder cancer has been demonstrated to have widespread genomic instability arising from defects in the DNA damage response and/or enhanced replication stress [57].These alterations facilitate the clonal evolution of tumour cells by accumulating driver aberrations, notably CNV and mutation [58].Molecular CNV profiles were notably heterogeneous between high-risk tumours (98 aberrations, comprising 60 amplifications covering 3795 genes and 38 deletions covering 5255 genes) and low-risk tumours (93 aberrations, comprising 49 amplifications covering 2965 genes and 44 deletions covering 7163 genes) (Figure 4A-D).Extensive heterogeneity in somatic  mutations was also detected in the two groups; Figure 4E shows the most prevalent mutational signatures in the samples.
In addition to the molecular features, the link between the signature and the pathways that are key drivers of bladder cancer was investigated.Strikingly, risk score presented on the rise with the increase of stromal activity pathway, e.g., pan-fibroblast TGFβ response signature (Pan-F-TBRS) and epithelial-mesenchymal transition (EMT) as well as DNA damage response pathways [59] (such as homologous recombination, mismatch repair, and nucleotide excision repair) (Figure 4F).Further investigation revealed that multiple tumourigenic pathways (cell cycle, extracellular matrix receptor interaction, focal adhesion, pathways in cancer, and Wnt signalling pathways) were remarkably enriched in high-risk tumours, with the notable enrichment of metabolism pathways (arachidonic acid metabolism, ascorbate and aldarate metabolism, drug metabolism cytochrome P450, linoleic acid metabolism, and porphyrin and chlorophyll metabolism) in low-risk tumours (Figure 4G and H), establishing the heterogeneity in pathway activity between low-and high-risk tumours.The signatures were computed as follows: where β i denotes the coefficient of LASSO regression for the ith gene and exp i denotes the expression value of the ith candidate gene.

Heterogeneity in antitumour immunity between low-and high-risk tumours
The immune cell infiltration was quantified using multiple analytical algorithms to ensure unbiased results.Interestingly, CD8+ T cells and activated natural killer (NK) cells were relatively more abundant in low-than in high-risk tumours (Figure 5A).In addition, a higher expression of common immune checkpoints (e.g.PDCD1 [60] and CD96 [61]) and immunostimulators (e.g.TNFRSF14, TNFSF13, and KLRK1) was detected in low-risk tumours (Figure 5B).These data demonstrated a difference in antitumour immunity between the two groups.

Distinct responses of low-and high-risk groups to ICB and commonly applied chemotherapeutic agents
In accordance with the lower HRD score and higher TMB and SNV neoantigens (Figure 6A-C), it was predicted that lowrisk patients would exhibit greater responses to ICB than high-risk patients.Interrupting one or more events within the cancer-immunity cycle facilitates escape from immunosurveillance.The RNA modification-related consensus signature-derived risk score displayed positive connections to several events, including the release of tumour antigens and the recruitment of basophils, eosinophils, neutrophils, and Th17 cells, decreased as the recruitment of Th2 cells and regulatory T (Treg) cells increased (Figure 6D), reflecting the associations of the signature with tumour responsiveness to ICB.The TIDE algorithm was also employed for ICB response estimation.Low-risk samples had higher dysfunction scores and lower exclusion and TIDE scores, indicating that this population was more likely to respond to ICB (Figure 6E-H).Together with the above-mentioned indicators, low-risk patients may benefit from ICB therapy.
Additionally, the IC 50 values of the commonly used chemotherapeutic agents were estimated.High-risk tumours presented remarkably lower IC 50 values for cisplatin and paclitaxel than low-risk tumours (Figure 6I and J), suggesting that high-risk patients are more likely to be suitable for chemotherapy.

Specific expression of genes from the RNA modification-related consensus signature in distinct T lymphocytes of bladder tumours
Next, we focused on the expression of genes from the RNA modification-related consensus signature in diverse intratumoural T cells isolated from bladder tumours.In total, 30,604 T lymphocytes from seven patients with bladder cancer were clustered into 18 cell clusters (Figure 7A).Six T cell types were marked in accordance with known marker genes: conventional CD4+ T (CD4Tconv), CD8+ T (CD8T), exhausted CD8+ T (CD8Tex), NK, proliferative T and Treg cells (Figure 7B and C).AKR1B1, ANXA1, CCNL2, OAS1, PTPN6, and TNFRSF14 were expressed in all T cell types, whereas SPINK1 was specifically expressed in CD4Tconv, CD8T, CD8Tex, and Treg cells (Figure 7D-J), indicating that these genes may mediate T cell-based antitumour immunity.
Transcriptional and post-transcriptional control of genes from the RNA modificationrelated consensus signature The transcriptional and post-transcriptional regulation of genes from the RNA modification-related consensus signature was further evaluated.Figure 8A illustrates transcriptional factors that potentially regulate the expression of AKR1B1, ANXA1, CCNL2, OAS1, PTPN6, SPINK1, and TNFRSF14.In addition, AKR1B1, ANXA1, CCNL2, OAS1, and PTPN6 were post-transcriptionally modulated by miRNAs (Figure 8B-D).

Discussion
To overcome the limitations of multigene signatures in underutilised data, inappropriate machine learning methods, and lack of rigorous validation in various cohorts and clinical trials, we adopted integrated multiple machine learning algorithms and their combinations to establish an RNA modification-related consensus signature [35,62].In this study, based on five independent cohorts, the best prediction signature was RSF with seven feature variables (PTPN6, AKR1B1, TNFRSF14, OAS1, ANXA1, SPINK1, and CCNL2), achieving the highest mean C-index of 0.656.This signature had the highest C-index in mean values of the four validation cohorts.Subsequently, according to the performance of the RSF-derived feature variables weighted by the corresponding regression coefficients acquired from a Cox model, patients with bladder cancer were classified into low-or high-risk groups according to the median risk score.Low-risk patients exhibited more favourable OS outcomes, ROC curves, and C-index than high-risk patients in the in-house and validated cohorts, demonstrating that the consensus signature had stable and reliable efficacy in prognostication.
Several clinical parameters affect prognostic outcomes in bladder cancer patients.While the incidence of bladder cancer is more significant in men, women have higher mortality than men, which results from disparate exposure to risk factors and the difference in importance of the sex steroid hormone pathway [63].AJCC staging has long served as a clinical guide in prognostic management and risk assessment, but inevitable limitations due to the heterogeneity of bladder cancer may lead to potential undertreatment or overtreatment [6].This study compared the efficacy of the consensus signature with that of conventional clinical parameters (age, sex, stage, and grade).The consensus signature presented notably superior accuracy in prognostic estimation according to the C-index assessment, further suggesting that it could act as a potential predictive marker for evaluating the prognosis of bladder cancer in clinical practice.
Over the last few decades, platinum-based chemotherapy has been regarded as the first-line treatment for eligible patients with bladder cancer, and the objective response rate of this therapeutic regimen is approximately 65-75 % [3].ICBs, such as programmed death/programmed death-ligand 1 (PD-1/PD-L1), cytotoxic T-lymphocyte-associated protein 4, and indoleamine-2, 3-dioxygenase-1, can promote the immune system to recognise and suppress the basic molecular targets of tumour cells [64].In the USA, PD-1/PD-L1 targeted drugs including atezolizumab and pembrolizumab, have been approved as first-line therapy for platinum-ineligible patients with bladder cancer and serve as maintenance therapy for those who are responsive to platinum [65].However, only a small proportion of patients with bladder cancer benefit from ICB therapy.In this study, indicators such as HRD, TMB, SNV, and TIDE were used to estimate ICB responsiveness.At the same time, IC 50 was used to evaluate the sensitivity to chemotherapy drugs.The results showed that low-risk patients might benefit from ICB, and high-risk patients are more likely to be suitable for chemotherapy, demonstrating that the consensus signature is an effective biomarker to predict patient responsiveness to medication and is helpful in guiding clinical therapeutic decision-making.
Based on multi-omics data, we further investigated the mutation and CNV features of the RNA modification-related consensus signature.Bladder cancer has been demonstrated to have widespread genomic instability arising from defects in the DNA damage response and/or enhanced replication stress [57].This study discovered that molecular CNV profiles and somatic mutations were notably heterogeneous between high-and low-risk tumours, which could account for their different prognoses.In addition to the molecular features, the link between the signature and the pathways that are key drivers of bladder cancer was investigated.The results showed that the risk score was positively correlated with stromal activity pathways.The dysregulation of the Pan-F-TBRS and EMT signalling systems is linked to the pathogenesis of cancers, such as cancer progression, chemoresistance, invasiveness, and metastasis [66,67].Further investigation revealed that multiple tumourigenic pathways were remarkably enriched in high-risk tumours, with a notable enrichment of metabolic pathways in low-risk tumours.These results may explain the poor prognosis in high-risk patients.
CD8+ T and NK cells are components of the lymphoid compartment in tumours [7].A previous study demonstrated that high CD8+ T cell infiltration generally predicts a better immunotherapy response [68].In this study, CD8+ T cells and activated NK cells were relatively more abundant in low-than in high-risk tumours.Additionally, a higher expression of common immune checkpoints and immunostimulators was detected in low-risk tumours.These results indicate that low-risk patients tend to have a better immune response, which is consistent with our conclusions.The results showed that the RNA modification-related consensus signatures, including AKR1B1, ANXA1, CCNL2, OAS1, PTPN6, SPINK1, and TNFRSF14, were specifically expressed in distinct T lymphocytes of bladder tumours at the single-cell level and were transcriptionally and post-transcriptionally modulated, indicating that these genes mediate T cell-based antitumour immunity and may be potential targets of immunotherapy.
This study has some limitations.First, our study was retrospective; further prospective studies and basic research are essential to supplement the relevant details of this study.Second, biased results are inevitable because the data analysis was based on public databases.Therefore, additional data from patients with bladder cancer are required to verify the practicality and accuracy of the consensus signature.Finally, the functions of most genes derived from the consensus signature remain unclear and need to be further explored through in vivo and in vitro experiments.

Figure 1 :
Figure 1: Multi-omics landscape of RNA modifiers in bladder cancer, and RNA modification-based consensus clusters.(A) Expression levels of RNA modifiers in bladder cancer and normal specimens.(B) Somatic mutation of RNA modifiers across bladder cancers.Mutation frequency is shown in the right panel, and mutation classification is marked by a unique colour.(C) Genomic location of RNA modifier copy number variations (CNVs).The numbers in the outermost circle represent chromosomes, while the innermost circle displays RNA modifiers.(D) CNV frequency of RNA modifiers across bladder cancers.Blue, copy number gain; green, copy number loss.(E) Correlation of RNA modifiers among samples.Red, positive association; blue, negative association. (F) Protein-protein interaction network based on RNA modifiers, and enrichment of Kyoto Encyclopedia of Genes and Genomes pathways.(G) Univariate Cox regression results of RNA modifiers with patient survival.Yellow, protective factor; blue, risk factor.(H) Ambiguously clustered pair (PAC) values at diverse optimal cluster number (k) values.The lowest PAC indicates the optimal number of clusters.(I) Tracking plot of the classification of samples under distinct k values.The smaller number of samples displaying different colours under distinct k values demonstrates relatively stable clusters.(J) Principal component analysis for the transcriptional difference of two clusters.(K) Consensus heatmap at k=2. White to blue indicates no clustering to always clustering.(L) Expression of prognostic RNA modifiers across the two clusters.(M) Survival probability of the two clusters.

Figure 2 :
Figure 2: Integrative establishment of an RNA modification-related consensus signature.(A) Sample dendrogram (upper) and heatmap of the two RNA modification clusters (lower).(B) Scale independence and mean connectivity values under a gradient of soft-threshold powers.(C) Gene dendrogram (upper) and co-expression modules based on expression similarity of genes (lower).(D) Relationships of co-expression modules with RNA modification classification.Correlation coefficient, and p-value are placed on box.(E) Relationships of module membership in the turquoise module with gene significance for RNA modification classification.(F) Multiple prediction signatures using the leave-one-out cross-validation framework and calculation of the concordance-index of each signature in the four validation cohorts.

Figure 3 :
Figure 3: The RNA modification-related consensus signature shows superiority to clinical variables in prognosis estimation.(A-E) Overall survival probability of low-and high-risk groups in the TCGA-BLCA, GSE32548, GSE48075, GSE69795, and GSE13507 cohorts.(F-H) One-, three-and five-year area under the receiver operating characteristic curve values in the above cohorts.(I) Concordance (C)-indexes of the signature in the above cohorts.(J, K) Comparison of C-indexes of the signature with those of traditional clinical parameters in the TCGA-BLCA, and GSE32548 cohorts.

Figure 4 :
Figure 4: Heterogeneity in genetic alterations, and signalling pathways between low-and high-risk tumours.(A, B) Copy number amplifications and deletions in high-risk tumours.The significance threshold was set as p<0.25 (marked by green line).(C, D) Copy number amplifications and deletions in low-risk tumours.(E) Somatic mutation landscape across low-and high-risk tumours.(F) Correlation analysis on the RNA modification-related consensus signature-derived risk score with well-established pathways.(G, H) Gene set enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes pathways in high-or low-risk tumours.

Figure 5 :
Figure 5: Heterogeneity in antitumour immunity between low-and high-risk tumours.(A) Immune cell infiltration landscape derived from multiple known algorithms both in high-and low-risk tumours.(B) Expression of well-defined immunomodulators in the two groups.

Figure 6 :
Figure 6: Distinct responses of low-and high-risk patients to immune checkpoint blockade and commonly used chemotherapeutic agents.(A-C) Comparison of homologous recombination deficiency (HRD) and tumour mutational burden (TMB) scores, and single nucleotide variant (SNV) neoantigens between low-and high-risk tumours.(D) Correlation analysis of the signature-derived risk score with cancer-immunity cycle activity.(E-G) Comparison of dysfunction, exclusion and tumour immune dysfunction and exclusion (TIDE) scores between groups.(H) Distribution of TIDE-predicted non-responders and responders across the two groups.(I, J) Comparison of estimated half-maximal inhibitory concentration (IC 50 ) values of cisplatin and paclitaxel between groups.*p<0.05;**p<0.01;***p<0.001.

Figure 7 :
Figure 7: Evaluation of the expression of genes from the RNA modification-related consensus signature in distinct T cells of blader tumours.(A) Clustering analysis of 30,604 T lymphocytes from seven patients with bladder cancer.(B) Identification of T cell types based on known marker genes.(C) Distribution of the proportions of each T cell type among seven bladder tumours.(D-J) UMAP plots illustrating the mapping of AKR1B1, ANXA1, CCNL2, OAS1, PTPN6, SPINK1 and TNFRSF14 in single T cells as well as their specific expression in diverse T cell types.

Figure 8 :
Figure 8: Transcriptional and post-transcriptional control of genes from the RNA modification-related consensus signature.(A) Transcriptional factorgene network.Circle represents genes from the RNA modification-related consensus signature, while rhombus represents transcriptional factors.(B-D) MicroRNA (miRNA)-gene networks.Circle indicates genes from the signature, while square indicates miRNAs.

Table  :
The comparison of clinical features between primary tumours and normal tissues.