Identification of m6A-related prognostic lncRNAs in KIRP
The 289 tumor samples and 32 normal tissues from KIRP patients were obtained and downloaded from the TCGA database and rearranged. After identifying 14086 lncRNAs in the TCGA-KIRP dataset, we extracted the expression profiles of 23 m6A-methylation regulators in the TCGA-KIRP dataset and then evaluated the relationship between the 23 m6A methylation regulators and 14086 lncRNAs. LncRNAs associated with one or more of the 23 m6A methylation regulators (|Pearson R| > 0.3 and p < 0.001) were defined as m6A-related lncRNAs, and the results showed that 1387 m6A-related lncRNAs were discerned as m6A-related lncRNAs. The co-expression networks of 23 m6A methylation regulators and lncRNAs were drawn using the Cytoscape software 3.7.0 and shown in Figure 1A. Univariate Cox regression analysis was performed to analyze m6A-related prognostic lncRNAs with a p-value of 0.001. The results showed that 32 m6A-related lncRNAs were associated with the overall prognosis of KIRP patients (Figure 1B and Table 1). Differences between tumor and normal tissues in 32 m6A-related prognostic lncRNAs were analyzed and displayed as a heatmap (Figure 1C). The correlation of 32 m6A methylation regulators and m6A-related lncRNAs were plotted using the Sankey diagram in Figure 1D. Furthermore, the expression levels of 32 m6A-related prognostic lncRNAs differed in KIRP tissues and normal tissues (Figure 1E).
Consensus clustering of m6A-related lncRNAs determined in KIRP
Consensus clustering was used to classify KIRP patients into subgroups based on the similarities shown by the expression of m6A-related lncRNAs. K= 2 showed the least overlap and the best cluster stability from k= 2 to 9, with the lowest CDF values; thus, lncRNAs were classified into cluster 1 and cluster 2 (Figures 2A-2C). To further evaluate the prognostic value of m6A-related lncRNAs based on the lncRNAs subtype survival analysis, the overall survival of KIRP patients in cluster 2 was lower than that in cluster 1 (p= 0.011, Figure 2D). Heatmap showed the differences in m6A-related prognostic lncRNAs expression in the two clusters and their relationship to clinicopathological parameters. M6A-related prognostic lncRNAs in cluster 2 were markedly associated with advanced pN stage, pM stage, pT stage, clinical stage, and grade in KIRP patients, and these lncRNAs were independent of age and gender (Figure 2E). Therefore, these consensus clustering results showed a significant correlation between m6A-related lncRNAs and the prognosis of KIRP patients.
Correlation analysis of m6A-related lncRNAs with the TME and immune cell infiltration
To investigate the role of m6A-related prognostic lncRNAs in the KIRP immune microenvironment, we first analyzed the differences in immuneScore and immune cell infiltration and exhibited a vioplot between cluster 1 and cluster 2 (Figure 3A). Immune cells such as memory B cells, resting CD4 memory T cells and resting mast cells were highly aggregated in cluster 1 (p< 0.05). In contrast, naive memory B cells, T cell CD8, activated CD4 memory T cells, follicular helper T cells, and macrophages M1 were highly clustered in cluster 2 (p< 0.05). The average immuneScore, stromalScore, and ESTIMATEScore in cluster 2 were significantly higher than in cluster 1 (Figures 3B-3D, all p< 0.01). Next, we examined the mRNA levels of immune checkpoints in each cluster and their association with m6A-related prognostic lncRNAs. The expression levels of PD-L1, HAVCR2, LAG3, PDCD1LG2, SIGLEC15, and TIGIT were prominently higher in cluster 2 than in cluster 1 (Figures 3E-3J). Furthermore, the expression of CTLA4, TIGIT, and SIGLEC15 was considerably elevated in KIRP tissues compared with normal tissues (Figure 3K). Interestingly, PD-L1, HAVCR2, LAG3, SIGLEC15, and TIGIT were negatively correlated with multiple m6A-related prognostic lncRNAs, and only PDCD1LG2 was positively correlated with the expression of some m6A-related lncRNAs (Figure 3L).
Construction and validation of the m6A-related prognostic lncRNAs signature
To further identify the most potent prognostic signatures, we performed a LASSO Cox regression analysis to screen and integrate 14 m6A-related prognostic lncRNAs to predict overall survival in KIRP patients (Figures 4A-4C). The 289 KIRP patients were randomly divided into training cohorts and test cohorts, and then the risk scores were calculated for each KIRP patient and divided into high-risk and low-risk groups. Survival curves revealed that the high-risk group of KIRP patients had a worse prognosis than those of the low-risk group in both the training cohort (Figure 4D, p < 0.001) and the test cohort (Figure 4G, p = 0.025). To assess the accuracy of the model in predicting survival in KIRP patients, the ROC results showed a curve (AUC) of 0.929, 0.921, and 0.930 for 1-, 2-, and 3-year overall survival rates, respectively, in the training cohort (Figure 4E). ROC curves also showed that m6A-related prognostic lncRNAs accurately predicted overall survival in the test cohort, with AUCs of 0.759, 0.766, and 0.754 for 1-, 2-, and 3-year overall survival rates, respectively (Figure 4H). The ROC curve and concordance index indicated that m6A-related prognostic lncRNAs were promising for predicting overall survival in KIRP patients. Subsequently, we achieved a risk curve and assessed the survival status and risk of m6A-related lncRNAs based on this curve. These features consisted of 14 m6A-related prognostic lncRNAs with excellent distinguishing performance in predicting the prognosis of KIRP patients. As risk scores increased, the number of deaths and the high-risk ratios enhanced. In addition, the expression of protective lncRNAs (ZKSCAN7-AS1, DSG2-AS1, AL133355.1, AC135050.6 and ADAMTS9-AS1) decreased with increasing risk score, while the expression of risk lncRNAs (CASC8, LUCAT1, AC099850.4, FALEC, AC096531.2, FAM66C, FOXD2-AS1, LINC01559, and AC010186.3) increased expression (Figures 4J-4K). These findings suggested that this m6A-related prognostic lncRNAs signature has robust and stable predictive power.
m6A-related lncRNA prognostic signature was an independent prognostic factor in KIRP
Univariate Cox regression analysis revealed that gender, clinical stage, and risk score were associated with the prognosis of KIRP patients in the training cohort (Figure 5A, all p< 0.01). Multivariate Cox regression analysis showed that the risk score was still significantly correlated with the prognosis of KIRP patients in the training cohort (Figure 5B, p = 0.003). In the test cohort, both univariate and multivariate regression analyses indicated that stage and risk scores were markedly correlated with prognostic factors in patients with KIRP (Figures5C-5D). Subsequently, based on these independent prognostic factors, a nomogram was developed to quantify the prediction of individual survival probability at 1, 2, and 3 years. The C-index for the nomogram associated with clinical parameters was 0. 868 in the training cohort (Figure 5E). The calibration curves showed that the predicted overall survival was consistent with the actual observations at three years in the training cohort (Figure 5G) and the test cohort (Figure 5H). Also, a nomogram of 1-, 3-, and 5-year overall survival rates based on m6A-related prognostic lncRNAs were determined using the rms package, and the C-index for the nomogram was 0. 888 in the training cohort (Figure 5F). Calibration curves showed agreement between predicted outcomes and actual survival, and similar results were obtained in the test cohort (Figures 5I-5J). The predictive power of the nomogram was similar to that of the ROC curve, indicating that the m6A-related prognostic lncRNAs model has high reliability. Furthermore, we performed the stratified survival analysis to assess the predictive power of m6A-related prognostic lncRNAs in KIRP patients with different clinicopathological parameters. In each subgroup, high-risk patients (especially men or women aged <= 60 years old) had lower overall survival than low-risk patients (Figures 6A-6C). Among KIRP patients with tumor grade 1-2 (Figure 6D, p= 0.045), tumor grade 3-4 (Figure 6E, p= 0.011), pT1-2 (Figure 6F, p= 0.012), pT3-4 (Figure 6G, p= 0.030), pN0 (Figure 6H, p= 0.007), pN1 (Figure 6I, p= 0.041), pM0 (Figure 6J, p< 0.001), stage I-II (Figure 6K, p= 0.045) and stage III-IV (Figure 6L, p= 0.011), the high-risk group had a worse prognosis than the low-risk group. These clinical data analyses demonstrated that m6A-related lncRNA prognostic signaturehad a good predictive performance.
m6A-related prognostic lncRNAs signatures correlated with clinicopathology and TME immune activity
To further disclose the correlation between clinical featuresand the prognostic risk models, the heatmap showed the high-risk and low-risk groups differed in pN stage, pM stage, pT stage, clinical stage, grade, gender, and cluster subtype (Figure 7A). The risk score also improved markedly as the clinical stage increased from stage I-II to stage III-IV (Figure 7B, p= 5.4e-09). Likewise, KIRP patients with pT3-T4, pN1-N2, and M1 had higher risk scores than T1-T2, N0, and M0 (Figures 7C-8E, all p< 0.01). The men, high-grade, and cluster 2 had higher risk scores than females and low-grade (Figures 7F-7G, all p< 0.01), and high-risk patients tended to be classified into cluster 2 (Figure 7H). These findings suggested that risk scores involving m6A-related prognostic lncRNA signatures were significantly associated with clinical characteristics in patients with KIRP.
We then comprehensively investigated the correlation between the risk score of m6A-related prognostic lncRNAs and immune cell infiltration. By combining differential and correlation analyses, we identified three types of immune cells correlated with risk scores involving m6A-related prognostic lncRNA signatures. The abundance of M2 macrophages was negatively correlated with the risk score (Figure 8A, p= 0.00059), and the abundance of M1 macrophages (Figure 8B, p= 0.00021) and plasma cells (Figure 8C, p= 0.049) was positively associated with the risk score. The expression level of immune checkpoint genes in high-risk and low-risk groups was also explored, the expression levels of HAVCR2 (p= 0.039) decreased, and LAG3 (p= 0.04) were increased in the high-risk group compared with the low-risk group (Figures 8D, 8E). These data indicated that the m6A-related lncRNAs prognostic features were associated with the KIRP immune microenvironment. Finally, we compared drug susceptibility to TKIs, chemotherapy, and targeted drugs, including sunitinib, axitinib, sorafenib, rapamycin, pazopanib, temsirolimus, and gemcitabine between risk groups. Risk stratification was significantly associated with sensitivity to sunitinib, rapamycin, pazopanib, and gemcitabine (Figures 8F-8L). These data showed that the risk scores of m6A-related lncRNAs prognostic signature was the potential targets of multiple drugs and had important significance for guiding the treatment of KIRP patients.