3.1 Genetic and transcriptional landscape of five types of RNA modification “writers” in BCa
In accordance with published articles depicting RNA modification [16, 23, 98, 99], a catalog of 34 RNA modification “writers” were enrolled into our study, including 9 m6A modification “writers” (METTL3, METTL14, RBM15, WTAP, KIAA1429, ZC3H13, METTL16, CBLL1 and RBM15B), 2 m6Am modification “writers” (PCIF1 and METTL4), 5 m1A modification “writers” (TRMT6, TRMT61A, TRMT10C, TRMT61B, RRP8), 15 APA modification enzymes (CPSF1/2/3/4, NUDT21, CPSF6/7, CSTF1/2/3, WDR33, FIP1L1, CLP1, PCF11, PABPN1), and 3 A-to-I modification enzymes (ADAR, ADARB1 and ADARB2) (Table S1).
To delineate genetic landscape of RNA modification “writers” in BLCA, we evaluated the frequency of non-silent somatic mutations in 34 “writers” based on mutational information of the TCGA-BLCA database. Specifically, 127 of 412 BLCA cases (30.83%) experienced mutations of RNA modification “writers”. METTL3 displayed the greatest mutation frequency (4%), followed by PCF11 (4%), KIAA1429 (3%), and WDR33 (3%). While the mutation frequency of ADARB1, METTL16 and CPSF7 was extremely low in BCa samples. Missense mutation constitutes the predominant type of mutations for each writer (Fig. 1a).
We demonstrated that BCa patients with mutant “writers” exhibited a significantly prolonged OS than those without mutation (HR = 0.53, 95% CI: 0.38–0.75, P < 0.001) (Fig. 1b), indicating that genetic alteration of “writers” potentially exerts a functional effect towards BCa tumorigenicity. GSEA was carried out to decipher biologic themes specific for “writers” mutated (“Writers” MUT) group (N = 126) and “writers” wild-type (“Writers” WT) group (N = 285) of patients in the TCGA-BLCA. “Writers” WT group was markedly enriched in carcinogenic activation pathways, such as angiogenesis, PI3K signaling, MAPK activity, P53 signaling, Jun kinase activity, and canonical Wnt signaling pathway (Fig. 1c; Table S2 and S3). Hence, the mutation of “writers” is prone to trigger functional alterations with prognostic significance in BCa.
We also investigated CNV alteration frequency of these “writers” and unraveled that ADAR, ADARB2, CLP1, and CPSF7 had a relatively high frequency of CNV amplification, while ZC3H13, RBM15B and RRP8 experienced a widespread frequency of CNV deletion (Fig. 1d). To determine whether CNV plays a considerable role in the expression of RNA modification “writers” in BCa patients, we attempted to assess the mRNA level of “writers” between normal and BCa tissues in the TCGA database. As depicted in Fig. 1e to 1i, a large proportion of enzyme-associated genes displayed relatively greater mRNA expression in BCa than normal tissues, highlighting the profound function of these “writers” in the occurrence and development of BCa. Moreover, RNA modification “writers” with CNV gain (such as ADAR, CLP1, and CPSF6) were significantly up-regulated in BCa samples than normal tissues. On the contrary, the expression of “writers” genes with CNV loss (including ZC3H13 and RRP8) was significantly diminished in BCa versus normal bladder tissues. Notably, certain “writers” (such as ADARB2 and PCF11) with widespread frequency of CNV gain harbored decreased mRNA level in BCa compared to adjacent non-tumor tissues (Fig. 1d, 1g and 1h).
To further elucidate the association between CNV values and mRNA expression in BCa samples, we stratified the TCGA-BLCA cohort into three groups according to CNV values of four “writers” characterized with exceeding 5% of CNV loss in BCa tissues, including CNV gain, CNV loss and non-significant alteration of CNV. Concretely, ZC3H13, RBM15B, RRP8, and RBM15 with CNV gain exhibited dramatically enhanced mRNA level than that with CNV loss, respectively. Nevertheless, the mRNA levels of above “writers” were significantly decreased or without remarkable alteration in CNV loss group compared with those in non-tumor samples (Figure S2a). Thus, CNV alteration partially explains why there is differential mRNA expression between tumor and normal samples [100]. Additional parameters, such as DNA methylation and transcription factors, are also endowed with the robust capacity to orchestrate gene expression in tumorigenesis [101, 102].
3.2 Prognosis and immune characteristics of RNA modification “writers” in BCa
Pairwise correlation analysis demonstrated that not only RNA modification “writers” in the same functional category exhibited a significant correlation in expression, but also a significant correlation was presented among mRNA levels of different category of “writers”. For example, BCa samples with high expression of A-to-I “writer” gene ADAR were accompanied by increased mRNA levels of eight m6A “writer” gene, including METTL14, RBM15, WTAP, VIRMA, ZC3H13, METTL16, CBLL1 and RBM15B, indicating a potential crosstalk between m6A and A-to-I modification in BCa (Fig. 2a). Whether co-expression phenomenon of these “writer” genes hints a functional correlation is a topic that motivates us to pursue further investigation. Additionally, prognosis analysis demonstrated that seven of 34 RNA modification “writers” were prognostic parameters of BCa cases in the TCGA-BLCA. Concretely, BCa patients with high METTL4 expression had a significantly longer survival time than those with low expression (HR = 0.71, 95% CI: 0.53–0.96, P = 0.0249), while those with higher KIAA1429 expression had a shorter survival time (HR = 1.35, 95% CI: 1.01–1.82, P = 0.0447) (Figure S2b). To further comprehensively expound the expression pattern of 34 “writers” in BCa, 1410 BCa patients from eight independent GSE sets were combined into a meta-GEO group in our study (Table S4).
RNA modification has recently gathered researchers’ attention because of its pertinent role in immune modulation and the TME. Thus, we reasoned the association between RNA modification “writers” and tumor immunity. As revealed in Fig. 2b and Figure S2c, METTL4 expression displayed a positive correlation with the subpopulation of tumor-infiltrating CD8+ T cells and natural killer cells (NK cells) with anti-tumor activity while was negatively associated with the proportion of myeloid derived suppressor cells (MDSCs) with immunosuppressive effect based on the meta-GEO cohort, and this is one of reasons why METTL4 were characterized with favorable prognostic value in BCa. Additionally, BCa patients with high KIAA1429 expression were characterized with an increased proportion of macrophage and Type 17 T helper cells (Th17 cells). We also mined the GSCALite web server and found that METTL4 might play a positive role in cell cycle and DNA damage response signaling pathways in BCa; nevertheless, there was no significant association between METTL4 expression and an array of tumor-promoting signals, including PI3K/Akt, TSC/mTOR, RAS/MAPK and receptor tyrosine kinase (RTK) pathways. Inversely, KIAA1429 expression was negatively related to apoptosis, and positively associated with PI3K/Akt and RTK pathways (Fig. 2c). The IHC experiment result also verified the greater expression of KIAA1429 protein in BCa tumor samples than normal tissues (Figure S2d). These results further implied the carcinogenic properties of KIAA1429 in BCa.
3.3 WGCNA used for the screening of METTL4
Considering that the broad impacts of m6Am “writer” METTL4 on epigenetic modifications and the positive association between METTL4 and multiple immune cells in BCa samples, we further emphatically discussed the prognostic significance of METTL4 in BCa. Initially, a total of 5657 prognosis-associated genes were extracted in 411 BCa patients, among which 3497 genes were associated with favorable prognosis (HR < 1, P < 0.05) and 2160 genes were related to unfavorable prognosis (HR > 1, P < 0.05). To select pivotal hub genes associated with BCa progression, above 5657 prognosis-associated genes were applied to cluster analysis by the “WGCNA” R package. On the basis of the standard scale-free network distribution, we carefully set the soft threshold power value as 7 to formulate a hierarchical clustering tree (dendrogram) of 5657 genes (Fig. 2d). According to the dynamic tree cut algorithm, the least gene number of each module and the minimum cut height was 50 and 0.25, respectively. The correlation of characteristic genes in integrated modules was over 0.9. We identified six co-expression modules containing all genes based on their degree of connectivity. The gray section represented background genes that did not belong to any modules (Fig. 2e). We ultimately assessed the correlation between MEs and clinical traits, including TNM stage, histological grade and survival status. Specifically, the brown module was characterized with the strongest negative correlation with survival status (r = -0.82, P < 0.0001), which was considered as the most significant module to the prognosis of BCa (Fig. 3a). There were four RNA modification “writers” in the brown module, of which MM value and the absolute value of GS of METTL4 ranked first (P < 0.0001) (Fig. 3b). Thus, METTL4 can be defined as one hub gene significantly related to survival status and BCa prognosis in the brown module.
GSEA analysis further revealed that six immune-associated pathways, such as regulation of receptor binding, Th17 type immune response, leukocyte adhesion to vascular endothelial cell, B cell differentiation, and negative regulation of IL-6 production and response to IFN-γ, were significantly enriched in high METTL4-expressing group (Fig. 3c and Table S5). Eventually, we also investigated whether METTL4 expression could predict patients’ response to ICB treatment. Specifically, 348 patients in IMvigor210 cohort displayed varying degrees of clinical benefits from anti-PD-L1 treatment, including complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). We demonstrated that patients with high METTL4 expression had a notably prolonged OS and prominently therapeutic advantage in IMvigor210 cohort (HR = 0.72, 95% CI = 0.54–0.96, P = 0.02) (Fig. 3d). Similar tendency was also identified in Snyder UC cohort (HR = 0.48, 95% CI = 0.18–0.79, P = 0.014) (Fig. 3e). The hypotheses regarding METTL4 expression as favorable prognostic indicator affected by its antitumor property came to light from the above bioinformatic approaches. Future mechanistic investigations of biological functions of METTL4 in BCa are required.
3.4 Immune landscape of RNA modification-associated patterns in BCa
On the basis of the expression profiles of 34 RNA modification “writers” in the meta-GEO cohort, we conducted unsupervised consensus clustering to stratify BCa patients with qualitatively varying RNA modification patterns into two distinct clusters, eventually including 1063 cases in Cluster 1 and 347 cases in Cluster 2, respectively (Fig. 4a and Figure S3a). Specifically, Cluster 1 had significantly increased presence of RBM15, WTAP, KIAA1429, TRMT6, TRMT61B, CPSF3, CPSF4, WDR33 and ADARB1, while Cluster 2 was characterized with elevated level of METTL14, ZC3H13, METTL16, PCIF1, METTL4, RRP8, CPSF1 and ADARB2 (Fig. 4b and Figure S3b). Furthermore, based on above-identified clusters, BCa patients in Cluster 1 and Cluster 2 were visibly separated into two discrete groups using three dimensional PCA, emphasizing that BCa cases are well stratified in line with the mRNA levels of 34 RNA modification “writers” (Fig. 4c). Survival analysis for two primary RNA modification subtypes demonstrated that compared with Cluster 1 modification pattern, Cluster 2 pattern was linked to significantly prolonged survival (HR = 0.63, 95% CI = 0.43–0.91, P = 0.013) (Fig. 4d). We further conducted GSVA to investigate the biological behaviors of above different RNA modification patterns. The carcinogenic activation pathways, including the epidermal growth factor activated receptor activity, TLR-2 signaling pathway, chemokine CXCL2 production and cell adhesion, were enriched relative to Cluster 1, indicating an inflammation activation and tumorigenesis status in Cluster 1. While Cluster 2 represented a metabolic or biosynthetic activation phenotype, prominently enriched in pathways related to the cyclic nucleotide catabolic process, CAMP-dependent protein kinase activity and DNA replication (Fig. 4e and Table S6).
We also unraveled the discrepancies concerning the compositions of tumor-infiltrating immune cells between two RNA modification clusters. Significant difference in immune cell fractions in two primary patterns were summarized in Table S7 and S8. As revealed in Fig. 4f and Figure S3c, Cluster 1 was characterized with an increased proportion of MDSCs with formidable immunosuppressive property (P < 0.0001) and Th17 cells (P < 0.0001). Conversely, cases in Cluster 2 exhibited prominent infiltration of activated CD8+ T cells with pronounced antitumor activity (P < 0.001), effector memory CD4+ T cells (P < 0.0001), and central memory CD8+ T cells (P = 0.032). Consistently, compared with cases in Cluster 2, those in Cluster 1 had significantly increased levels of MDSC marker genes (including STAT2, S100A9, CXCL2, CSF1, PTGS2, TREM1, CEBPB) while significant down-regulation of activated CD8+ T cell marker genes (such as CD8A, IFN-γ, IL-13, and FASLG) (Fig. 4g). In summary, RNA modification patterns exert an effect on the proportions of infiltration by specific immune cell types while fail to change the types of infiltrating immune cells.
3.5 Establishment of the RMS model and its clinical significance in BCa
We further determined 632 RNA modification phenotype-associated DEGs between Cluster 1 and Cluster 2. The biological processes with significant enrichment associated with these DEGs were enriched in purine nucleotide metabolic process, methylation, DNA replication, and cell cycle G1/S phase transition, all of which were closely related to RNA processing (Table S9 and Figure S4a). In consideration of the heterogeneity and complexity of RNA modification, we attempted to establish a risk score system named the RMS (RNA Modification “Writers” Score). Firstly, we confirmed that 119 of 632 DEGs was significantly correlated with OS through univariate Cox regression analysis (Table S10). To reveal potential DEGs with the optimal prognostic performance, we utilized LASSO Cox analysis, and 14 DEGs were incorporated into our subsequent analysis. Notably, IFNLR1 (HR = 1.95, 95% CI: 1.47–2.58, P < 0.0001) was identified as a crucial RNA phenotype-associated and prognosis-related DEGs with the largest beta value (0.206), implying its optimal prognostic efficiency in BCa patients (Table S11; Figure S4b and S4c). Furthermore, we performed normalization of the expression levels of 14 DEGs in the TCGA and meta-GEO cohort with average value = 0 and SD = 1, thus acquiring a uniform cutoff value as stratified standard. Then, we quantified the RNA modification status of each BCa patient by weighting the normalized mRNA level of each DEG to the regression coefficient. The concrete formula was the following: RMS = 0.2062 * Expr IFNLR1 + 0.1822 *ExprPCDHB11 + (-0.1428) * Expr TIMM21 +……+ 0.0057 * Expr CRELD1+0.0028* Expr FOXG1. Ultimately, we calculated the RMS for each BCa case in the meta-GEO and stratified all cases into RMS-high and -low cohorts based on the median value (3.344) (Table S12).
As showed in Fig. 5a, the RMS classified BCa cases with high or low risk score into two discrete sections, highlighting that the RMS distribution of BCa cases in the low-risk group was greatly different from those with high risk score. There was a high degree of consistency among the risk score distribution, the heatmap of 14 prognostic DEGs’ expression and survival status of BCa case in the meta-GEO cohort (Fig. 5b). Notably, the cutoff point (3.344) also served as a classification indicator in the TCGA-BLCA cohort. Kaplan-Meier curve revealed that high RMS was significantly correlated with more unfavorable clinical outcome of BCa cases in the meta-GEO (HR = 3.00, 95% CI: 2.19–4.12, P = 1.06e− 11) (Fig. 5c) and TCGA-BLCA cohort (HR = 1.53, 95% CI: 1.13–2.09, P = 0.006) (Fig. 5d). Additionally, patients with high RMS had a shorter disease-specific survival (DSS) than their RMS-low counterparts in GSE32894 (HR = 18.3, 95% CI: 4.30–77.6, P < 0.001) (Figure S4d). In GSE31684 dataset with recurrence data, the RMS was significantly negatively associated with recurrence-free survival (RFS) (HR = 2.02, 95% CI: 1.06–3.88, P = 0.033) (Figure S4e). We also investigated the correlation between the RMS and cluster classifier to evaluate the RMS model’s accuracy. As revealed in Fig. 5e, the RMS of BCa samples in Cluster 1 was significantly higher than that of cases in Cluster 2 (P = 0.025). We found that 579 out of 742 (78.03%) samples with high RMS were overlapped with the samples in Cluster 1, and 184 out of 668 (27.54%) cases in RMS-low group overlapped with the individuals in Cluster 2 (Figure S4f and S4g). Therefore, there is a high degree of coincidence among three computational methods of classification.
Univariate Cox analysis indicated that certain clinical variables, including age, M classification, N classification, T classification, TNM stage, and RMS exhibited an impact on the survival of BCa patients (Fig. 5f). Above significant parameters were included into subsequent multivariate Cox regression analysis. The corresponding findings revealed that age > 65 years (HR = 1.87, 95% CI: 1.33–2.64, P = 0.00032), advanced N classification (HR = 1.62, 95% CI: 1.17–2.25, P = 0.00403), and high RMS (HR = 1.61, 95% CI: 1.17–2.24, P = 0.00388) remained adverse and independent prognostic factor in BCa (Fig. 5f). These findings imply that the RMS model is independent of conventional clinical variables and can predict the survival of BCa with comparatively satisfactory performance. Additionally, based on WGCNA method, the green module was most significantly positively correlated with survival status of BCa (r = 0.76, P < 0.0001) (Fig. 3A), and IFNLR1 was one of the most significant hub genes in the green module (Fig. 5g), indicating that IFNLR1 was one of the DEGs that was the most significantly correlated with prognosis than other 13 DEGs in the RMS model. Additionally, we performed IHC of BCa samples and matched bladder tissues from 84 cases with primary BCa to further investigate IFNLR1 protein expression level in BCa. IFNLR1 protein was significantly increased in BCa tissues compared with adjacent normal bladder samples (Fig. 5h).
To confer physicians with a visualized approach to predict the long-term survival of BCa patients, the nomogram model encompassing the RMS signature and significant clinical risk factors identified by multivariate Cox analysis was formulated. As illustrated in Fig. 6a, N classification made the greatest contributions to risk points, followed by age and the RMS model. The C-index of the nomogram was 0.745 (95% CI: 0.697–0.798) under 1000 bootstrap replication. The calibration curves for the OS probability of 1-, 3-, and 5-year in BCa cases demonstrated a good agreement between nomogram prediction and practical observation (Fig. 6b). The time-ROC curves were established to compare the predictive efficiency of this nomogram with that of N classification, age and the RMS. For the ROC curve of 1-year survival, the AUC of nomogram (0.671) was higher than that of age (0.622), N classification (0.582), and the RMS (0.54) (Fig. 6c). The nomogram to predict 3-year OS obtained the optimal AUC of 0.675, followed by N classification (0.603), age (0.588), and the RMS (0.558) (Fig. 6d). The AUC for the nomogram, N classification, age and the RMS to predict 5-year survival were 0.706, 0.625, 0.584, and 0.601, respectively (Fig. 6e). The decision curve analysis (DCA) of the nomogram was characterized with the optimal net benefits, followed by N classification, age, and the RMS (Fig. 6f). In sum, the nomogram incorporating N classification, age and the RMS exhibits a relatively satisfactory predictive performance for the long-term survival of BCa patients.
3.6 Molecular subtypes and functional annotation associated with the RMS in BCa
We further illustrated the functional characteristics of the RMS signature through analyzing the association between the RMS model and known biological processes-associated gene sets identified by MSigDB [87], emphasizing that high RMS was significantly associated with stromal activation status and cancer progression-associated pathways, such as inflammatory response, NF-KB-mediated TNF-a, epithelial-mesenchymal transition (EMT), angiogenesis, IL-6/JAK/STAT3 signaling (Fig. 7a and Figure S4h).
Based on a comprehensive molecular subtypes landscape of UC established by Sjödahl et al.’s study [70], UC cases were classified into five molecular subtypes, including urobasal A (MS1 subdivided into MS1a and MS1b), genomically unstable (MS2a subdivided into MS2a1 and MS2a2), urobasal B (MS2b2.1), squamous cell carcinoma (SCC)-like (MS2b2.2), and one highly infiltrated by nontumor cells (MS2b1). Notably, above molecular subtypes differed in survival patterns in which urobasal A exhibited favorable prognosis, genomically unstable and the infiltrated group were with moderate prognosis, and the urobasal B and the SCC-like were characterized with the shortest survival [70]. We compared the difference of the RMS among above five molecular subtypes through analyzing data downloaded from GSE32894. The SCC-like subtype showed the highest RMS, followed by the urobasal B, genomically unstable, the infiltrated subgroup and urobasal A (Fig. 7b). Additionally, there was a significant discrepancy in the distribution of molecular subtypes between RMS-high and -low group. The urobasal A subtypes was primarily clustered in RMS-low group, conversely, the SCC-like and the infiltrated subtype were markedly concentrated on RMS-high group (Fig. 7c).
To further decipher potential biological processes associated with different molecular subtypes of BCa patients in GSE32894, GSVA was performed implying that tumorigenesis- associated biological processes were significantly enriched in the SCC-like and the urobasal B subtype, including WNT, NOTCH, angiogenesis, and IL2-STAT5 signaling pathways. In contrast, the biological pathways activated in the urobasal A subtype were significantly correlated with the heme metabolism, protein secretion and peroxisome (Fig. 7d). Consistently, the SCC-like and the urobasal B-related signaling pathways were prevailingly enriched in RMS-high cases, while the enrichment score of urobasal A-related biological processes were markedly clustered in RMS-low group (Figure S4h). BCa patients in the urobasal B and the SCC-like subtype were prone to be diagnosed at more advanced stage compared with those in the urobasal A subtype (Fig. 7e), which was also significantly correlated with diminished OS (HR = 12.3, 95% CI: 1.36–111, P = 0.026 for the urobasal B subtype) (Fig. 7f). Previous results in our report demonstrated that the RMS positively correlated with BCa patients’ degree of malignancy. Thus, high RMS roughly corresponding to the urobasal B and the SCC-like subtype indicates unsatisfactory prognosis, which is potentially partly ascribed to the activation of EMT, WNT, angiogenesis, and additional signaling pathways mediating BCa tumorigenicity and tumor metastasis.
3.7 Pan-cancer analysis of RMS model-associated genes
Initially, we explored the correlation between CNV and mRNA expression in 14 RMS model-associated genes in 33 kinds of tumors and revealed that CHMP7 expression was significantly modulated by CNV in almost all cancers, followed by SEPHS1 and AASDHPPT (Figure S5a). Specifically, a majority of RMS model-associated genes was characterized with heterozygous amplification of CNV in adenoid cystic carcinoma, while homozygous amplification was prone to occur in OV, esophageal cancer, and unconditioned stimulus (Figure S5b). Thus, the findings highlight that the CNV of RMS model-associated genes is various among different tumors and it is essential to investigate the source of the heterogeneity. Furthermore, we explored the difference in mRNA expression between tumor and normal sample and revealed that the fold difference in the expression of RMS model-associated genes was the greatest in LUSC. Concretely, GDPD5 and IL28RA were significantly downregulated in BRCA than normal samples, while ROMO1 was overexpressed in BRCA samples (Figure S5c). Additionally, pathway analysis demonstrated that RMS model-associated genes principally triggered cell cycle and EMT pathway while exerts an inhibitory effect on apoptosis and RAS/MAPK pathway (Figure S5d). Therefore, our RMS model genes potentially plays a crucial role in malignant progression of tumors.
3.8 Difference in post-transcriptional events between RMS-high and -low groups in BCa
To elucidate the functional effect of RNA modification “writers” on post-transcriptional characteristics of BCa patients, we investigated APA events of each gene in the TCGA-BLCA. Initially, we analyzed APA alterations between 246 BCa cases with high or low RMS and determined the prognostic significance of transcripts with significant 3’UTR alterations. A total of 11598 APA events remained for differential analysis, and there were 503 genes with significantly lengthened 3’UTR (ΔPDUI > 0.1) and 96 transcripts with markedly shortening 3’UTR (ΔPDUI < 0.1) in RMS-high group, respectively (P < 0.05) (Fig. 8a and Table S13), and shortening APA events in RMS-high group were characteristic with worse OS based on univariate Cox regression analysis (Fig. 8b and Table S14), thus indicating that usage of a PAS may exacerbate BLCA malignancy. Specifically, the transcripts of CCNO (ΔPDUI = -0.16, P = 0.003) and PAOX (ΔPDUI = -0.15, P = 0.03) both underwent statistically significant 3’UTR shortening in patients with high RMS, which was associated with worse survival in BLCA (HR = 1.92, 95% CI: 1.30–2.86, P = 0.001 for CCNO; HR = 1.52, 95% CI: 1.02–2.22, P = 0.039 for PAOX) (Fig. 8c). A report have demonstrated that CCNO is overexpressed in cervical squamous cell carcinoma (CSCC) and RACK1/ miR-302b/c/d-3p-mediated CCNO inhibition can dampen the progression of CSCC [103]. Suppression of PAOX is sufficient to widen the therapeutic index of cytotoxic drugs and overwhelm DNp73-mediated chemoresistance in cancers [104]. Thus, we speculate that shortening 3’UTR of CCNO and PAOX in BLCA samples with high RMS potentially results in loss of several RNA regulatory elements, such as miRNA binding sites, thus enabling the upregulation of oncogenes expression and the progression of BLCA.
3.9 Potential role of the RMS in antitumor chemotherapy and antibody-drug conjugates (ADC) therapy
To further assess the relationship between the RMS and drug response of BCa cell lines, we determined 34 significantly correlated pairs between the RMS and drug response in the GDSC database based on Spearman correlation analysis [95]. Specifically, there was significant correlation between drug sensitivity and the RMS in 8 pairs, including EGFR inhibitor HG-5-88-01 (Rs = -0.804, P = 0.005), CSF1R inhibitor GW-2580 (Rs = -0.43, P = 0.016), and AR inhibitor Bicalutamide (Rs = -0.383, P < 0.0001). Conversely, 26 pairs displayed drug resistance significantly related to the RMS, including JNK1 inhibitor ZG-10 (Rs = 0.867, P < 0.0001) and CDK9 inhibitor THZ-2-49 (Rs = 0.625, P < 0.0001) (Fig. 9a and Table S15). Additionally, we also explored the potential signaling pathways of drug-targeted genes. As revealed in Fig. 9b, drugs whose sensitivity was linked to high RMS primarily targeted hormone-related, ADCK4, and EGFR signaling pathways, while drugs whose resistance was related to high RMS mostly targeted DNA replication, cell cycle and PI3K/MTOR signaling pathways. Thus, above findings indicate that RNA modification patterns are related to drug response of tumors. The RMS potentially develops into a novel biomarker to confer a reference for appropriate clinically interventional strategies.
ADCs are novel targeted agents that concatenate a cytotoxic drug (also known as cytotoxic payload or warhead) by a linker to a monoclonal antibody (mAb) which can specifically reach target antigens expressed on cancer cellular surface and deliver a potent cytotoxic payload to the tumor location, ultimately strengthening the chemotherapeutic efficacy and minimizing toxicity to normal tissue. The target antigen should be abundantly expressed on tumor cells while is not expressed or at a low level in normal tissues in an ideal setting, thus lowering off-target toxicity [105]. Currently, certain ADCs have been approved by the US Food and Drug Administration (FDA) for the cancer therapy (Table S16) [106–113]. Two target antigens (ERBB2 and TROP2) were lineage-specific markers of two out of above approved ADCs——Trastuzumab deruxtecan, and Sacituzumab govitecan, which have consistently high expression across the BCa tumor population than normal samples in the TCGA-BLCA (Fig. 9c). We further evaluated the differences in the expression of seven target antigen molecules of ADC in RMS-low and -high groups. The target antigens, including ERBB2 and TROP2, were preferential expression on RMS-high BCa samples with a relative low expression on RMS-low subgroup (Fig. 9c). Together, above findings implied that RNA modification patterns are potentially associated with ADC sensitivity.
3.10 Predictive value of the RMS in immunotherapeutic efficacy
Immunotherapies of blocking T-cell inhibitory molecules PD-L1 and PD-1 have undoubtedly emerged a significant breakthrough in anticancer intervention. Meanwhile, it is urgent for us to make judgment about which subset of patients can benefit most from immunotherapies [114]. Therefore, we investigated the predictive power of the RMS for patients’ response to ICB therapy based on two immunotherapy cohorts. Patients with low RMS exhibited significantly prolonged OS than those with high RMS in IMvigor210 cohort (HR = 0.76, 95% CI: 0.58–0.99, P = 0.040) (Fig. 10a). For IMvigor210 cohort, Chi-squared test demonstrated that compared with RMS-high group, RMS-low group was endowed with markedly increased proportion of the sum of CR and PR patients while significantly diminished the sum of PD and SD cases (P = 0.037) (Fig. 10b). Likewise, CR patients were characterized with the lowest RMS compared with their counterparts with other types of response (Fig. 10c). Significant therapeutic advantage and strengthened clinical response to anti-MAGE-A3 immunotherapy in patients with low RMS were also confirmed in Montoya melanoma cohort (Fig. 10d and 10e). Additionally, we also validated the predictive performance of the RMS in anti-MAGE-A3 immunotherapy, with a satisfactory AUC value of 0.712 (Fig. 10f). Collectively, cases with lower RMS are more possibly to reap better prognosis and enhanced clinical benefit from ICB therapy.
Accumulated evidence has emphasized that patients with elevated TMB, higher neoantigen burden, certain DNA repair mutations, mismatch repair deficiency, and higher PD-L1 expression level are correlated with improved objective response, durable clinical benefit, and prolonged long-term survival when receiving ICB therapy [115, 116]. Based on tumor-associated immune phenotypes depicted in IMvigor210 cohort, patients with low RMS were characterized with significantly increased PD-L1 level (Fig. 11a). Similarly, cases in RMS-low group had significantly strengthened TMB and neoantigen burden than those with high RMS (Fig. 11b and 11c), indicating a potential response to ICB. Patients with the combination of low RMS and high TMB/neoantigen burden displayed the optimal survival advantage (HR = 0.51, 95% CI: 0.33–0.79, P = 0.003 for Low RMS with high TMB; HR = 0.48, 95% CI: 0.31–0.76, P = 0.002 for Low RMS with high neoantigen burden) (Fig. 11d and 11e). We further explored the difference in the RMS among three phenotypes, including “immune inflamed”, “immune excluded”, and “immune desert” [117]. As illustrated in Fig. 11f, patients with an immune-inflamed phenotype exhibited the lowest RMS compared with the other two phenotypes. Above findings partly explain why immunotherapy is prone to exert intensive antitumor effect in the low RMS subset. Our aforementioned results also demonstrated that MDSC which is recognized to mediate immune tolerance in the TME was significantly activated in RMS-high group, indicating that RMS-high tumors potentially represent “cold tumors” with resistance to immunotherapy. Furthermore, AUC value evaluating the capacity of the RMS model, TMB, TNB and PD-L1 to differentiate responders from non-responders was 0.677 (95% CI = 0.589–0.765), 0.652 (95% CI = 0.549–0.755), 0.690 (95% CI = 0.595–0.785), and 0.625 (95% CI = 0.517–0.733), respectively. The results also illustrated that the RMS in combination with TMB, TNB and PD-L1 had the optimal predictive power, with the highest AUC of 0.828 (95% CI = 0.714–0.941), followed by TMB combined with TNB and PD-L1 (AUC = 0.797, 95% CI = 0.678–0.916), the RMS combined with TNB (AUC = 0.765, 95% CI = 0.671–0.859), the RMS combined with TMB (AUC = 0.742, 95% CI = 0.641–0.843), and the RMS combined with PD-L1 (AUC = 0.708, 95% CI = 0.595–0.822) (Fig. 11g). Briefly, these results may introduce the novel piece to the atlas of RNA modification patterns’ influence on the efficacy of immunotherapy.