Identification of prognostic and m6A-related lncRNAs in BLCA
A total of 14,086 lncRNAs were screened from TCGA data, and the detailed clinicopathological information of patients is shown in Additional file 1: Table S1. We determined that 745 m6A-related lncRNAs were significantly correlated with 23 m6A-related genes using Pearson correlation analysis (Fig. 1a). After excluding patients without cancer or survival data, we merged the survival data with the lncRNA expression data of individual patients (final patient number = 403). Subsequently, we identified 51 lncRNAs related to OS (Additional file 2: Table S2). Figure 1b shows that the expression of these prognostic m6A-related lncRNAs differed significantly between normal and BLCA tissues.
Consensus clustering of m6A-related lncRNAs identified in molecular subtypes of BLCA
Based on the expression profile of the prognostic m6A-related lncRNAs, we categorized the patients with BLCA into two groups: cluster 1 (n = 130) and cluster 2 (n = 273) (Fig. 2a–c). Survival analysis demonstrated that patients in cluster 2 had a worse OS than those in cluster 1 (p = 0.010, Fig. 2d). Figure 2e shows the relationships between the two clusters, revealing a distinct difference between the clusters in terms of the tumor stage (p < 0.05). Thus, the results of the consensus clustering analysis are associated with the progression of BLCA and survival of patients with BLCA.
GSEA and immune-related analysis of the two clusters
The results of GSEA show that terms associated with multiple tumor hallmarks, such as apoptosis, epithelial-mesenchymal transition, inflammatory response, TNFA signaling, and IL6-JAK-STAT3 signaling, were predominantly enriched in cluster 2 (NOM p and FDR q-value < 0.001, Fig. 3a). C2 KEGG analysis similarly revealed several tumor-related signaling pathways and cell adhesion activities in cluster 2, whereas only three gene sets were significantly enriched in cluster 1, with metabolic pathways at a NOM p-value < 0.05 (Fig. 3b–e). Surprisingly, several enriched immune-related signaling pathways and genes were identified through both Hallmark and C2 KEGG analyses. C7 collection analysis indicated that multiple immune functional gene sets were enriched in clusters 1 and 2 (Fig. 3f).
Interestingly, we found a significant difference between the scores in clusters 1 and 2 (Fig. 4a). In addition, among the TICs in the BLCA, regulatory T cells (Tregs), neutrophils, and M2 macrophages were significantly associated with the clusters (Fig. 4b–c). These findings suggest that the expression of m6A-related lncRNAs in the clusters might have immunomodulatory effects on the TME. Given this, we attempted to determine the correlation between the clusters and some immune checkpoints. It is worth noting that seven types of immune checkpoint molecules, except for GAL9, were highly expressed in cluster 2 (Fig. 4d). The correlations between the 51 lncRNAs and the eight immune checkpoint molecules are shown in Additional file 3: Fig. S1 and most of these correlations were significant. Based on these results, we speculate that the poor prognosis in cluster 2 was probably due to the upregulation of the immune checkpoint molecules.
Construction and validation of the m6A-RLPS
There were no significant differences between the training cohort(n = 203) and validation cohort (n = 200) in terms of any of the clinicopathological parameters (Additional file 4: Table S3). Next, LASSO Cox regression analysis incorporated nine prognostic m6A-related lncRNAs to the m6A-RLPS for predicting the OS of patients with BLCA (Fig. 5a–c). The risk score for each patient was calculated as follows: risk score = (− 0.433 × AC020911.1 expression) + (0.059 × KCNQ1OT1 expression) + (− 0.359 × AC104532.2) + (− 0.841 × AC006160.1) + (− 0.239 × EHMT2-AS1 expression) + (0.242 × AC097359.2 expression) + (0.768 × AP001469.1 expression) + (− 0.039 × AC007686.3 expression) + (− 0.048 × AL022322.1). Univariate Cox regression analysis supported the remarkable prognostic significance of all nine lncRNAs, among which AC097359.2, AP001469.1, and KCNQ1OT1 were identified as risk factors with a hazard ratio (HR) > 1, whereas the others were found to be protective factors (all p values < 0.01, Fig. 5d).
Next, we collected the expression data from the GEO DataSets platform for validation analysis of the 9 lncRNAs. Ultimately, only KCNQ1OT1 among the m6A-RLPS was detected in the GSE31189, GSE31684, and GSE51493 cohorts for subsequent analyses. The results showed that KCNQ1OT1 was markedly overexpressed in BLCA samples compared with normal tissues or in the high pathological T stage. Furthermore, Kaplan-Meier log-rank test supported the KCNQ1OT1 was a remarkable prognostic risk factor for BLCA in GSE31684 (Additional file 5: Fig. S2). These findings were in accordance with the results of TCGA analysis.
Subsequently, Kaplan–Meier curves showed that in the training cohort, patients in the low-risk group had an improved OS compared with those in the high-risk group (p < 0.001, Fig. 6a). Similar results were obtained in the validation cohort and the entire cohort (both p < 0.001, Fig. 6d, Additional file 6: Fig. S3a). The distributions of risk score, survival status, and corresponding lncRNA expression in each cohort are shown in Fig. 6b, Fig. 6e, and Additional file 6: Fig. S3b, respectively. The ROC curves showed that the m6A-RLPS could accurately forecast the OS in the training cohort, and the AUCs were 0.729, 0.707, and 0.769 for the 1-, 3-, and 5-year OS rates, respectively (Fig. 6c). Similar results were subsequently verified in the validation cohort and the entire cohort (Fig. 6f, Additional file 6: Fig. S3c). These findings indicate that the prognostic signature has a robust and stable predictive efficiency.
The m6A-RLPS was an independent prognostic indicator in BLCA
The result of the univariate analysis indicated that age (p < 0.01), tumor stage, and risk score (both p < 0.001) were closely associated with prognosis (Fig. 7a). Multivariate analysis validated that age, tumor stage, and risk score were independent prognostic factors in patients with BLCA (Fig. 7b). Stratification survival analysis showed that high-risk patients had an observably worse OS than low-risk patients in every subgroup (Fig. 7c–l), indicating a notable predictive performance of the m6A-RLPS. Nomograms for the 3- and 5- year OS rates based on the independent predictors determined from the multivariate analysis are shown in Fig. 7m. A certain point was generated for each covariate, and a total nomogram score, which was correlated with the 3- and 5-year OS rates, was calculated for every patient. The nomogram showed favorable accuracy in predicting the OS, with a C-index of 0.71 (95% CI: 0.61–0.77) and 0.68 (95% CI: 0.62–0.74) for the training and validation cohorts, respectively. Moreover, the calibration curves revealed that there was an appreciable agreement between the predictive outcome and actual survival, and similar conclusions were obtained in the validation cohort (Additional file 7: Fig. S4). In addition, we further determined the utility of the m6A-RLPS risk score across 33 kinds of TCGA cancer using the transcriptome and clinicopathological data acquired from the UCSC Xena project (http://xena.ucsc.edu). The univariate Cox regression results showed that the m6A-RLPS risk score was significantly related to OS in 9 types of cancer (BLCA, COAD, KICH, KIRC, LGG, PAAD, SKCM, STAD, and UVM). Among them, the m6A-RLPS risk score was a risk factor in BLCA, PAAD, SKCM, and STAD (HR > 1, p < 0.05), whereas all others were protective factors. Considering the possibility of death from non-tumor causes during follow-up, we also analyzed the relationship between the m6A-RLPS risk score and Disease-specific survival (DSS) in 33 TCGA tumors. The univariate Cox regression results also showed that the m6A-RLPS risk score was a risk factor in BLCA, PAAD, OV, SKCM, and UCSE, but was a protective factor in KIRC, LGG, and UVM (Additional file 8: Table S4). These results demonstrated the predictive value of the m6A-RLPS for some other cancers.
The m6A-RLPS was correlated with clinicopathology and TME immune activity
Figure 8a shows that the expression of lncRNAs included in the m6A-RLPS was significantly correlated with cluster type, immune score, tumor grade, tumor stage, T stage, and N stage (all p values < 0.01). In addition, the Student’s t-test demonstrated that the risk score increased with increasing immune score, clinical stage, T stage, and N stage (Fig. 8d–g) and that high-risk patient tended to be gathered in cluster 2 (Fig. 8h), but the risk score was not significantly correlated with age and sex (Fig. 8b-c). These findings suggest that the m6A-RLPS can influence the progression of BLCA.
By combining the difference and correlation analyses, we found that six types of TICs, including plasma cells, Tregs, M0 macrophages, M2 macrophages, activated dendritic cells, and neutrophils were strikingly associated with the m6A-RLPS (Fig. 9). Of these associations, three showed positive correlations (M0 macrophages, M2 macrophages, and neutrophils), whereas the others showed negative correlations. In addition, we discovered that the m6A-RLPS also had a strong positive correlation with the TME scores obtained using the ESTIMATE algorithm (Fig. 10a). We also noticed that the expressions of immune checkpoints, except for GAL9, were increased in high-risk patients and positively correlated with the risk score, reflecting the effect of the immune checkpoints on TME and poor oncological outcomes (Fig. 10b–c).
