RNA sequencing profiles (FPKM value) and matched clinical information of ccRCC, pRCC and chRCC were retrieved from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov) via the Genomic Data Commons (https://portal.gdc.cancer.gov/) utilizing TCGAbiolinks package [18]. The raw count value was converted to TPM value. Three RCC datasets were integrated and batch effects were removed through Combat function of sva package [19]. Finally, 128 normal samples and 893 RCC samples (chRCC (N = 65), ccRCC (N = 539) and pRCC (N = 289) were included in this study. Copy number alterations (CNAs) of RCC that were preprocessed by GISTIC algorithm were obtained from the UCSC Xena data portal (http://xena.ucsc.edu/). Additionally, RNA expression profiling of 101 ccRCC patients from the E-MTAB-1980 dataset was adopted for external validation. In total, 23 m${}^6$A regulators containing 8 writers (Cbl proto-oncogene like 1 (CBLL1), Vir like m${}^6$A methyltransferase associated (VIRMA), methyltransferase 14/3 (METTL14/3), RNA binding motif protein 15 (RBM15), RNA binding motif protein 15B (RBM15B), WT1 associated protein (WTAP), zinc finger CCCH-type containing 13 (ZC3H13)), 2 erasers (ALKBH5 and FTO) and 13 readers (ELAV like RNA binding protein 1 (ELAVL1), FMRP translational regulator 1 (FMR1), heterogeneous nuclear ribonucleoprotein A2/B1 (HNRNPA2B1), heterogeneous nuclear ribonucleoprotein C (HNRNPC), insulin like growth factor 2 mRNA binding protein 1/2/3 (IGF2BP1/2/3), leucine rich pentatricopeptide repeat containing (LRPPRC), YTH domain containing 1/2 (YTHDC1/2), YTH m${}^6$A RNA binding protein 1/2/3 (YTHDF1/2/3)) were collected from the published literature. The mRNA expression of the m${}^6$A regulators was compared between 128 normal and 893 RCC samples. Univariate-cox regression analyses were presented to investigate the correlation between the m${}^6$A regulators and RCC prognosis in the TCGA dataset. The results were visualized via forestplot package (version 2.0.1, https://cran.r-project.org/web/packages/forestplot) [20]. The expression of prognostic m${}^6$A regulators was validated utilizing the Human Protein Altas (https://www.proteinatlas.org/).

The m${}^6$A regulators were inputted into the STRING online database (http://string-db.org/) [21]. The interactions were visualized into a protein-protein interaction (PPI) network via Cytoscape software (version 3.8.2, https://cytoscape.org/) [22]. Pearson correlation analysis was performed to evaluate the correlation of the mRNA expression of the m${}^6$A regulators across RCC samples.

Through Tumor Immune Dysfunction and Exclusion (TIDE) approach, the response to immunotherapy was predicted following the tumor immune evasion mechanisms [21]. The expression similarity between subtypes and the patients who might respond to anti-PD-1 and anti-CTLA4 therapy was detected utilizing Subclass Mapping (SubMap) algorithm [23].

Unsupervised clustering analyses were applied to estimate the number of unsupervised classes across 893 RCC samples based on the expression profile of 23 m${}^6$A regulators by ConsensuClusterPlus package (version 1.60.0, https://bioconductor.org/packages/release/bioc/html/ConsensusClusterPlus.html) [24]. The classification accuracy was assessed via t-distributed stochastic neighbor embedding (t-SNE) method. Survival differences among clusters were compared by Kaplan-Meier curves and log-rank tests.

In total, 50 hallmark pathways were retrieved from the Molecular Signatures Database, which may comprehensively reveal the major biological functions of humans [25]. The enrichment levels of above pathways were quantified with GSVA package [26]. From the REACTOME dataset (https://reactome.org/) [27] and the Gene set enrichment analysis (GSEA) website (http://www.gsea-msigdb.org/gsea/index.jsp) [28], 32 oxidative stress-related genes were collected. The activity of oxidative stress was estimated with single-sample gene set enrichment analysis (ssGSEA) function [29].

Two stemness indicators: gene expression-based stemness index (mRNAsi) and DNA methylation-based stemness index (mDNAsi) were calculated in RCC samples by building the stemness index models with one-class logistic regression (OCLR) machine-learning method [30]. The stemness value was ranging from 0 (no gene expression) to 1 (complete gene expression).

Cancer immunity cycle contains release of cancer cell antigens (step 1), cancer antigen presentation (step 2), priming and activation (step 3), trafficking of immune cells to tumors (step 4), infiltration of immune cells into tumors (step 5), recognition of cancer cells by T cells (step 6), and killing of cancer cells (step 7) [31]. The activation of each step was quantified in RCC samples by ssGSEA function.

Immune score and stromal score were obtained using estimation of stromal and immune cells in malignant tumours using expression data (ESTIMATE) algorithm based on gene expression data, which represented the fractions of immune and stromal cells in RCC samples [32]. The enrichment levels of immune cells were estimated in RCC samples through the ssGSEA algorithm according to the expression matrix of gene symbols of tumor-infiltrating immune cells (TIICs) [33]. The mRNA expression of immune checkpoints was calculated in each RCC specimen.

Differentially expressed genes (DEGs) were screened between m${}^6$A clusters rough limma package [34]. Genes with adjusted p 0.05 were considered as DEGs. DEGs shared by any two clusters were identified as m${}^6$A-related DEGs.

Functional enrichment analyses of m${}^6$A-related DEGs were presented with clusterProfiler package (version 4.9.0.2, https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html), containing gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) [35]. GO categories comprised biological processes (BPs), cellular components (CCs) and molecular functions (MFs). Terms with adjusted p 0.05 were significantly enriched by m${}^6$A-related DEGs.

RCC patients were clustered into distinct gene clusters through applying unsupervised clustering analyses based on the extracted m${}^6$A-related DEGs. By using ConsensuClusterPlus package, the number of gene clusters and their stability were determined across RCC samples. The accuracy of the gene clusters was validated with t-SNE method. Survival differences among gene clusters were observed with Kaplan-Meier curves and log-rank tests. Univariate-cox regression analyses were carried out for extracting prognostic m${}^6$A-related DEGs with p 0.05 across RCC patients. These extracted DEGs were used for feature selection through recursive feature elimination (RFE) with random forest, followed by 10-fold cross-verification utilizing caret package. The m${}^6$A scoring system was quantified in individual tumors based on the curated expression profiles of the finally identified DEGs by conducting a principal component analysis (PCA) utilizing the Boruta algorithm, named the m${}^6$A score. Both principal component 1 and 2 were chosen as signature scores. The m${}^6$A score was quantified according to the following formula: m${}^6$A score = $\sum$ (PC1${}_{\text{i}}$ + PC2${}_{\text{i}}$), where i indicated the expression of m${}^6$A-related DEGs.

The mRNA expression profiles and follow-up data of patients with anti-PD-1 therapy were downloaded from the GSE78220 dataset (N = 27) in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) repository [36] and the literature reported by Liu et al. [37] (N = 121). The prognostic value of m${}^6$A score was externally validated in the two immunotherapy cohorts.

Pathways were probed by GSEA package [28]. The “c5.bp.v6.2.symbols.gm” gene set was set as the reference set, which was obtained from the molecular signatures database. Pathways with adjusted p 0.05 were significantly associated with m${}^6$A score.

By employing pRRophetic algorithm [38], the half-maximal inhibitory concentration (IC50) values of chemotherapy agents were estimated for RCCs on the basis of the genomics of drug sensitivity in cancer (GDSC; www.cancerrxgene.org/) cell line expression spectrum [39] and mRNA expression profiling.

Drug sensitivity data of human cancer cell lines (CCLs) were retrieved from the the Cancer Therapeutics Response Portal (CTRP) (https://portals.broadinstitute.org/ctrp) and Preservice Research Institute for Science and Mathematics (PRISM) database (https://depmap.org/portal/prism/) datasets. Both two datasets offer the area under the dose–response curve as an evaluation indicator for the responses to hundreds of compounds. Reduced area under the curve (AUC) indicates higher sensitivity to a specific drug. By K-nearest neighbor method, the missing AUCs were imputed. Since the CCLs in both datasets were retrieved from the cancer cell line encyclopedia (CCLE) project (https://portals.broadinstitute.org/ccle/) [40], expression profiling data in CCLE project were employed for CTRP and PRISM analyses.

Statistical analyses were conducted with R packages (version 4.0.2). Kaplan-Meier curves of overall survival (OS), disease free survival (DFS), disease-specific survival (DSS), progression free survival (PFS) and progression free interval (PFI) were depicted for RCC patients and survival differences between groups were compared with log-rank tests. Comparison between two groups was presented via student’s t test or Wilcoxon test. One-way ANOVA or Kruskal-Wallis test was used for conducting comparison between three or more groups. Two-sided p 0.05 indicated statistical significance.

This study collected three types of RCC ccRCC, pRCC and chRCC from TCGA datasets and removed batch effects after integration (Fig. 1A,B). Finally, 893 RCC samples (chRCC (N = 65), ccRCC (N = 539) and pRCC (N = 289) and 128 normal samples were included for further analyses. In total, 23 m${}^6$A regulators were collected and the prevalence of CNAs among these m${}^6$A regulators was analyzed in RCC samples. In Fig. 1C, YTHDC2 had the prevalent gain and IGFBP2, YTHDF2 and RBM15B had the widespread loss across RCC tissues. The mRNA expression of these m${}^6$A regulators was compared in RCC and normal specimens. We found that most m${}^6$A regulators were abnormally expressed in RCC compared to normal tissues (Fig. 1D), including writers (CBLL1, METTL14, ZC3H13), eraser (FTO) and readers (FMR1, HNRNPA2B1, HNRNPC, IGF2BP1/2/3, LRPPRC, YTHDC1/2, YTHDF2/3). We also compared the expression of m${}^6$A regulators in stage I&II and III&IV RCC tissues. As depicted in Fig. 1E, writers (CBLL1, METTL14, METTL3, RBM15B, ZC3H13), and readers (FMR1, IGF2BP1, IGF2BP3, LRPPRC, YTHDC1, YTHDF1) displayed remarkable differences between stage I&II and III&IV. Additionally, we noted that there were significant differences in the expression of writers (RBM15, WTAP), erasers (ALKBH5, FTO) and readers (HNRNPA2B1, HNRNPC, IGF2BP2/3, YTHDC1, YTHDF2) among three RCC types ccRCC, pRCC and chRCC (Supplementary Fig. 1). Through Submap algorithm, we predicted the responses to anti-PD-1 and anti-CTLA4 therapy. As a result, among three RCC types, ccRCC patients were more likely to respond to anti-PD-1 therapy. There was no significant difference in the response to anti-CTLA4 therapy among three RCC types (Supplementary Fig. 1). By univariate-cox regression analyses, prognostic implication of each m${}^6$A regulator was evaluated across RCC patients. METTL14, ZC3H13, FTO, LRPPRC and YTHDC1 were protective factors of RCC, while HNRNPA2B1 and IGF2BP1/2/3 acted as risk factors of RCC (Fig. 1F). Further multivariate cox regression analysis showed that METTL14 and LRPPRC were protective factors, while IGF2BP3 and HNRNPA2B1 were identified as risk factors, which was consistent with the results of univariate-cox regression analysis (Supplementary Fig. 2). Additionally, we evaluated the relationships between m${}^6$A regulators and DFS, DSS, and PFS outcomes. Our results demonstrated that the m${}^6$A regulators IGF2BP1/3 were significantly associated with DFS, DSS, and PFS, indicating that they might contribute to RCC progression (Fig. 1G–I). The PPI network revealed the tight interactions among the m${}^6$A regulators (Fig. 1J). Also, there were significant correlations between the m${}^6$A regulators at the mRNA levels across RCC samples (Fig. 1K). Above data highlighted the important implications of the interactions of m${}^6$A regulators in RCC progression.

Fig. 1.

Systematic analyses of genetic variation, expression, prognostic implication and interactions of N6-methyladenosine (m6A) regulators in renal cell carcinoma (RCC). (A,B) Before and after removing batch effects by integrating three types of RCC: clear cell RCC (ccRCC), papillary RCC (pRCC) and chromophobe RCC (chRCC) from the Cancer Genome Atlas (TCGA) datasets. PCA, Principal component analysis. (C) The frequency of copy number variation (CNV) (red: gain and blue: loss) of m${}^6$A regulators across RCC samples. (D) An overview of the mRNA expression of m${}^6$A regulators in RCC and normal specimens. p values were calculated by unpaired student’s t test, *p 0.05 and ****p 0.0001. Red: up-regulation and green: down-regulation. (E) The mRNA expression of m${}^6$A regulators in stage I&II and stage III&IV RCC cases. *p 0.05, ***p 0.001 and ****p 0.0001. (F–I) Univariate-cox regression analyses of the correlation between m${}^6$A regulators and RCC patients’ overall survival (OS), disease free survival (DFS), disease-specific survival (DSS), and progression free survival (PFS). (J) The protein-protein interaction (PPI) network of writers (green), erasers (blue) and readers (red). (K) The Pearson correlations between m${}^6$A regulators at the mRNA levels across RCC samples. Red: positive correlation and green: negative correlation, *p 0.05, **p 0.01 and ***p 0.001.

Using the human protein altas database, we verified the expression of prognostic m${}^6$A regulators in RCC and normal kidneys. As depicted in Fig. 2A, METTL14, ZC3H13, HNRNPA2B1 were lowly expressed in RCC than normal kidneys. Meanwhile, FTO, LRPPRC, YTHDC1, IGF2BP1/2/3 were highly expressed in RCC compared with normal kidneys. We also evaluated the interactions of m${}^6$A regulators with 32 oxidative stress-related genes in RCC. In Fig. 2B, m${}^6$A regulators was significantly linked to most oxidative stress-related genes, indicating that the activity of oxidative stress might be modulated by m${}^6$A modification during RCC progression.

Fig. 2.

Protein expression of prognostic m${}^{\mathrm{6}}$A regulators and interactions of m${}^{\mathrm{6}}$A regulators with oxidative stress in RCC. (A) Validation of the expression of prognostic m${}^6$A regulators in RCC and normal kidneys from the Human Protein Altas database. Scale bar, 200 µm. (B) Heatmap of the interactions of m${}^6$A regulators and oxidative stress-related genes in RCC. Red: positive correlation and green: negative correlation, *p 0.05; **p 0.01.

By performing consensus clustering analysis, we classified 893 RCC samples into five m${}^6$A modification subtypes based on the expression profile of the m${}^6$A regulators (Fig. 3A), named m${}^6$A modification subtype A (N = 85), B (N = 160), C (N = 105), D (N = 237) and E (N = 292). The t-SNE confirmed the accuracy of this classification (Fig. 3B). Survival analyses were carried out among the five m${}^6$A modification subtypes. In Fig. 3C, we found that there were significant differences in prognosis among subtypes, where subtype A exhibited the worst survival outcomes, subtype E and B had the prominent survival advantage, and subtype C and D had the moderate survival time. The mRNA expression of the m${}^6$A regulators was compared among subtypes. As a result, subtype C had the lowest expression of the m${}^6$A regulators, followed by subtype D (Fig. 3D). Most m${}^6$A regulators were up-regulated in subtype A, B and E. To validate the reliability of this m${}^6$A modification classification, the E-MTAB-1980 dataset was adopted. Consistently, ccRCC patients were clearly classified into five m${}^6$A modification subtypes (Supplementary Fig. 3). The notable survival difference among the subtypes was also proven (Supplementary Fig. 3). The activation of the major biological processes was quantified in each RCC specimen vis ssGSEA. Almost all pathways were inactivated in subtype C and D (Fig. 3E). However, most pathways were activated in subtype A, B and E. Cancer stemness was quantified by two indices: mDNAsi and mRNAsi. In Fig. 3F, subtype B had the lowest mDNAsi value while subtype E had the highest mDNAsi value. But mRNAsi value was the lowest in subtype E (Fig. 3G). The somatic copy number alteration (SCNA) was observed in each RCC specimen. Among subtypes, subtype A exhibited the highest SCNA level (Fig. 3H).

Fig. 3.

Construction of the m${}^{\mathrm{6}}$A modification subtypes with different prognosis, activation of pathways, cancer stemness and somatic copy number alteration (SCNA) across RCC samples. (A) Consensus clustering analyses for identifying the best clustering number (k = 5) across RCC samples based on the expression matrix of m${}^6$A regulators. (B) The t-SNE for assessing the clustering accuracy across RCC samples. (C) Survival analyses of RCC patients in different m${}^6$A modification subtypes with log-rank tests. (D) An overview of the mRNA expression of m${}^6$A regulators in each m${}^6$A modification subtype. Red: up-regulation and green: down-regulation. (E) An overview of the activation levels of the 50 main signaling pathways in RCC specimens from five m${}^6$A modification subtypes. (F) The mDNAsi score, (G) mRNAsi score and (H) SCNA level among different m${}^6$A modification subtypes. p values were determined with Kruskal-Wallis test.

By ssGSEA method, we quantified the infiltration levels of immune cells in RCC tissues. As a result, subtype A had the highest infiltration levels of most tumor-infiltrating immune cells, but subtype B displayed the lowest infiltration levels of most immune cells (Fig. 4A,B). With the same approach, the activity of oxidative stress was quantified. The significant heterogeneity in oxidative stress was found among five subtypes (Fig. 4C). Among them, subtype B had the lowest activity of oxidative stress. The mRNA expression of immune checkpoints PD-1 and PD-L1 was compared among subtypes. We found that patients in subtype A exhibited the highest mRNA expression of PD-1 (Fig. 4D) and PD-L1 (Fig. 4E). The enrichment levels of each step in the cancer immunity cycle were evaluated across RCC samples. Our results showed that most steps exhibited the highest activation in subtype A, while most steps exhibited the lowest activation in subtype B (Fig. 4F). The overall infiltrations of immune cells and stromal cells were calculated by ESTIMATE method. As expected, the highest immune score and stromal score were found in subtype A (Fig. 4G,H). Meanwhile, we found that there were the lowest immune score and stromal score in subtype B. Furthermore, subtype A displayed the lowest tumor purity while the highest tumor purity was detected in subtype B (Fig. 4I). Microsatellite instability (MSI) was quantified and compared among m${}^6$A modification subtypes. In Fig. 4J, our data showed that there was the lowest MSI level in subtype E but there was the highest MSI level in subtype B.

Fig. 4.

Assessment of the m${}^{\mathrm{6}}$A modification subtypes with distinct oxidative stress and tumor immunity across RCC samples. (A,B) The infiltration levels of immune cells among five m${}^6$A modification subtypes. (C) Activity of oxidative stress among five m${}^6$A modification subtypes. (D,E) The mRNA expression of immune checkpoints programmed cell death 1 (PD-1) and programmed cell death ligand 1 (PD-L1) among five m${}^6$A modification subtypes. (F) The activation of each step in the cancer immunity cycle among five m${}^6$A modification subtypes. (G–I) Immune score, stromal score, and tumor purity among five m${}^6$A modification subtypes. (J) The microsatellite instability (MSI) level among five m${}^6$A modification subtypes. p values were calculated with Kruskal-Wallis test, ns: not significant; *p 0.05; **p 0.01; ***p 0.001.

By comparing the DEGs between m${}^6$A modification subtypes, we finally identified 733 m${}^6$A-related DEGs (Fig. 5A; Supplementary Table 1). GO enrichment analyses indicated that these m${}^6$A-related DEGs were significantly correlated to methylation modification processes such as histone modification, covalent chromatin modification and histone H3-H4 methylation (Fig. 5B). Meanwhile, we found that the m${}^6$A-related DEGs were markedly enriched in mRNA modification pathways (such as mRNA surveillance pathway) and carcinogenic pathways (such as renal cell carcinoma, TGF-beta signaling pathway and central carbon metabolism in cancer; Fig. 5C), which confirmed that m${}^6$A modification played a non-negligible role in the tumor progression. Among 733 m${}^6$A-related DEGs, 512 were significantly associated with RCC prognosis according to univariate-cox regression analyses (Supplementary Table 2).

Fig. 5.

Establishment of the m${}^{\mathrm{6}}$A genomic subtypes with different survival outcomes across RCC samples. (A) Venn diagram for the m${}^6$A-related differentially expressed genes (DEGs) by comparing the DEGs between m${}^6$A modification subtypes. (B,C) Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment results of the m${}^6$A-related DEGs. The length of the rectangle indicated the number of enriched genes in a specific term. The more the color tended to red, the smaller the adjusted p value. (D) The four m${}^6$A genomic subtypes across RCC samples based on the m${}^6$A-related DEGs via unsupervised clustering analyses. (E) The t-SNE for the transcriptome profiles among the four m${}^6$A genomic subtypes. (F) Kaplan-Meier survival curves for the four m${}^6$A genomic subtypes across RCC samples. p values were determined with log-rank tests.

To further validated the regulation mechanism of m${}^6$A modification, we conducted unsupervised clustering analyses based on the obtained five m${}^6$A cluster-related genes. On the basis of the expression profiling of the 512 prognostic m${}^6$A-related DEGs, 893 RCC samples were clustered into four m${}^6$A genomic subtypes by unsupervised clustering analyses, named m${}^6$A genomic subtype A (N = 226), B (N = 318), C (N = 140) and D (N = 195; Fig. 5D). The t-SNE plots confirmed that there was a distinct difference on transcriptome profiles among the four m${}^6$A genomic subtypes (Fig. 5E). Prognostic analyses showed that there was a significant difference in survival outcomes among m${}^6$A genomic subtypes (Fig. 5F), where m${}^6$A genomic subtype B had the most prominent survival advantage and m${}^6$A genomic subtype A exhibited the poorest survival outcomes.

To explore the biological processes of the four m${}^6$A genomic subtypes, we used GSVA to analyze the activity of the top 50 signaling pathways. In Fig. 6A, m${}^6$A genomic subtype B and D displayed the activation of most pathways, while m${}^6$A genomic subtype A and C had the inactivation of most pathways. Among four m${}^6$A genomic subtypes, subtype D displayed the highest activity of oxidative stress (Fig. 6B). The roles of the m${}^6$A genomic subtypes in regulating the TME immune infiltration were further analyzed. We found that several antitumor immune cells such as activated CD8 T cells, activated CD4 T cells, central memory CD4 T cells, central memory CD8 T cells, effector memory CD4 T cells, effector memory CD8 T cells, activated dendritic cells and natural killer T cells displayed the relatively high infiltration levels in m${}^6$A genomic subtype A (Fig. 6C). Meanwhile, the most protumor immune cells such as immature dendritic cells, and plasmacytoid dendritic cells exhibited the low infiltration levels in m${}^6$A genomic subtype A and C as well as the high infiltration levels in m${}^6$A genomic subtype B and D. The activity of cancer immunity cycle was also analyzed. As a result, priming and activation of the immune system, recognition of cancer cells by T cells and killing of cancer cells had the highest activity in m${}^6$A genomic subtype A (Fig. 6D). Above steps exhibited the relatively low activity in m${}^6$A genomic subtype B and C. This indicated that m${}^6$A genomic subtype A might be an inflammatory phenotype, while m${}^6$A genomic subtype B and C might be non-inflammatory phenotypes. Also, we quantified the mRNA expression of immune checkpoints. As expected, PD-1 expression was the highest in m${}^6$A genomic subtype A, while its expression was the lowest in m${}^6$A genomic subtype C (Fig. 6E). Furthermore, m${}^6$A genomic subtype B had the highest PD-L1 expression, while m${}^6$A genomic subtype C had the lowest PD-L1 expression (Fig. 6F). These data indicated that m${}^6$A genomic subtype A and B might be more sensitive to anti-PD-1/PD-L1 therapy.

Fig. 6.

The m${}^{\mathrm{6}}$A genomic subtypes characterized by oxidative stress, tumor microenvironment (TME) cell infiltration and transcriptome features across RCC samples. (A) An overview of the activity of the major 50 signaling pathways in RCC specimens from the four m${}^6$A genomic subtypes. Red: activation and green: inactivation. (B) Activity of oxidative stress in the four m${}^6$A genomic subtypes. (C) The enrichment scores of tumor-infiltrating immune cells in the four m${}^6$A genomic subtypes. (D) The activity of each step in the cancer immunity cycle in the four m${}^6$A genomic subtypes. (E,F) The mRNA expression of immune checkpoints PD-1 and PD-L1 in the four m${}^6$A genomic subtypes. p values were calculated with Kruskal-Wallis test, ns: not significant; *p 0.05; **p 0.01; ***p 0.001.

Due to the heterogeneity and complexity of m${}^6$A methylation modification, a m${}^6$A scoring system was generated for quantifying m${}^6$A modification subtypes of each RCC patient based on the expression profiling of the prognostic m${}^6$A-related DEGs that were identified by PCA method using the Boruta algorithm, named m${}^6$A score (Fig. 7A). For better illustrating the clinical value of m${}^6$A score, this study analyzed the distribution of m${}^6$A score in different clinicopathological characteristics (Fig. 7B). Our data showed that patients with age $\ge$65 had the significantly lower m${}^6$A score compared to those with age 65. There was a significantly decreased m${}^6$A score in male patients than female patients. With the increase of grade, stage and T, m${}^6$A score was gradually reduced across RCC specimens. Metastatic patients had significantly lower m${}^6$A score than non-metastatic patients. Patients with N1-2 displayed markedly reduced m${}^6$A score in comparison to those with N0. In addition, we also explored the relationship between clinical factors and m${}^6$A score. Low m${}^6$A score was corresponding to more dead cases (p 0.001), advance stage and grade cases (Supplementary Fig. 4). Based on the median m${}^6$A score, RCC patients in the TCGA cohort were separated into high and low m${}^6$A score group. Survival analysis were then presented for comparing the prognostic divergence between groups. As a result, patients with high m${}^6$A score displayed the prominent advantage in OS, DSS and PFI in comparison to those with low m${}^6$A score (Fig. 7C–E). Consistently, better OS outcomes were observed in high m${}^6$A score group relative to low m${}^6$A score group in the E-MTAB-1980 dataset (Fig. 7F). The prognostic value of m${}^6$A score was also externally verified in the two anti-PD-1/PD-L1 immunotherapy cohorts. Both in the GSE78220 dataset and the dataset reported by Liu et al. [37], high m${}^6$A score indicated undesirable outcomes than low m${}^6$A score among patients who received anti-PD-1/PD-L1 therapy (Fig. 7G,H). Finally, we estimated the correlation between m${}^6$A cluster, genomic cluster and m${}^6$A score using sankey plot (Supplementary Fig. 5). Most of cases from m${}^6$A cluster E and genomic cluster B were attribute to high m${}^6$A score, which can partly explain why m${}^6$A cluster E and genomic cluster B associated with a better survival outcome.

Fig. 7.

Generation of m${}^{\mathrm{6}}$A scoring system and assessment and external validation of the prognostic value of m${}^{\mathrm{6}}$A score for RCC patients. (A) Identification of the final m${}^6$A-related DEGs that were used for construction of m${}^6$A score by principal component analysis (PCA) method utilizing Boruta algorithm. (B) The distribution of m${}^6$A score in different clinicopathological characteristics across RCC patients, including age, sex, stage, grade and T, N, M. p values were determined with Wilcoxon or Kruskal-Wallis tests. (C–E) Kaplan-Meier curves of OS, DSS and PFI for RCC patients with high and low m${}^6$A score in the TCGA cohort. (F) Kaplan-Meier curves of OS between high and low m${}^6$A score ccRCC patients in the E-MTAB-1980 cohort. (G,H) Kaplan-Meier curves of OS between high and low m${}^6$A score groups for patients treated with anti-PD-1/PD-L1 therapy in the GSE78220 dataset and the dataset reported by Liu et al [37].

The activity of the major signaling pathways was compared in high- and low-m${}^6$A score RCC patients. We found that most pathways such as carcinogenic pathways and inflammation-related pathways were distinctly activated in high-m${}^6$A score group compared to low-m${}^6$A score group. Meanwhile, in comparison to high-m${}^6$A score group, KRAS signaling, coagulation and myogenesis were distinctly activated in low-m${}^6$A score group (Fig. 8A). Consistently, GSEA results revealed that apoptosis, focal adhesion, insulin signaling pathway, MAPK signaling pathway, oocyte meiosis, pathway in cancer and WNT signaling pathway were significantly up-regulated in high-m${}^6$A score samples than low-m${}^6$A score samples (Fig. 8B). No pathways were significantly enriched in low-m${}^6$A score specimens. Furthermore, we evaluated the differences in tumor immunity between high- and low-m${}^6$A score groups. Most immune checkpoints such as LAIR1, HAVCR2 and ICOS displayed the higher mRNA expression in high-m${}^6$A score group compared to low-m${}^6$A score group (Fig. 8C). TNFRSF18, TNFSF14, TMIGD2, LAG3 and CD70 were significantly up-regulated in low-m${}^6$A score group than high-m${}^6$A score group. There were increased infiltration levels of central memory CD8 T cell, eosinophil, effector memory CD4 T cell, memory B cell, type 2 T helper cell, mast cell, neutrophil, natural killer cell, plasmacytoid dendritic cell and immature dendritic cell in high-m${}^6$A score group than low-m${}^6$A score group (Fig. 8D). Meanwhile, activated CD8 T cell, CD56dim natural killer cell, monocyte, MDSC and activated B cell exhibited the higher infiltration levels in low-m${}^6$A score group in comparison to high-m${}^6$A score group. Most oxidative stress-related genes had higher expression levels in high- than low-m${}^6$A score group (Fig. 8E).

Fig. 8.

Correlation between m${}^{\mathrm{6}}$A score and signaling pathways, tumor immunity and oxidative stress in RCC samples. (A) An overview of the correlation between m${}^6$A score and the enrichment score of the major 50 signaling pathways in RCC samples. The m${}^6$A scores from left to right gradually increased. Red indicated high activation of a specific pathway and green indicated low activation of a specific pathway. (B) Gene set enrichment analysis (GSEA) for the activated signaling pathways in high-m${}^6$A score samples. (C) An overview of the correlation between m${}^6$A score and the mRNA expression of immune checkpoints in the RCC cohort. Red represented high expression and green represented low expression for a specific immune checkpoint. (D) An overview of the correlation between m${}^6$A score and the infiltration levels of immune cells in RCC specimens. Red indicated high infiltration level and green represented low infiltration level for a specific immune cell. (E) The mRNA expression of oxidative stress-related genes in high and low m${}^6$A score groups. Comparisons between two groups were analyzed with Wilcoxon test, *p 0.05; **p 0.01; ***p 0.001.

IC50 values of four approved targeted drugs including axitinib, pazopanib, sorafenib and sunitinib were estimated in each RCC sample. The differences in IC50 values were compared between high- and low-m${}^6$A score RCC patients. In comparison to low-m${}^6$A score group, there were significantly reduced IC50 values of axitinib, pazopanib and sorafenib in high-m${}^6$A score group (Fig. 9A–C). This indicated that RCC patients with high-m${}^6$A score were more sensitive to axitinib, pazopanib and sorafenib. Conversely, lower IC50 value of sunitinib was found in low-m${}^6$A score group compared to high-m${}^6$A score group (Fig. 9D), demonstrating that low-m${}^6$A score patients were more likely to benefit from sunitinib. Potential therapeutic agents were predicted for low-m${}^6$A score patients. Two drug dataset (CTRP and PRISM) were adopted to identifying candidate CTRP and PRISM-derived agents with higher drug sensitivity in low-m${}^6$A score patients. Differential drug response analyses between high- and low-m${}^6$A score groups were carried out for predicting drugs with reduced AUCs in low-m${}^6$A score patients. As a result, six CTRP-derived compounds (CR-1-31B, daporinad, leptomycin B, methotrexate, oligomycin A and ouabain; Fig. 9E) and ten PRISM-derived compounds (bosutinib, gemcitabine, GSK1070916, halobetasol-propionate, mesna, OTX015, PD-0325901, Ro-4987655, romidepsin and vincristine; Fig. 9F) had reduced AUCs in low-m${}^6$A score group compared to high-m${}^6$A score group, indicating that they could become potential therapeutic agents against low-m${}^6$A score patients. Additionally, potential druggable targets were predicted for low-m${}^6$A score patients. As a result, we found that m${}^6$A score was positively associated with CERES scores of four druggable targets (NMU, MOCOS, UCHL1, and RARRES1), with negative relationships with their protein expression (Fig. 9G,H). Therefore, NMU, MOCOS, UCHL1, and RARRES1 might become potential druggable targets for low-m${}^6$A score patients.

Fig. 9.

Comparison of the response to drugs in high- and low-m${}^{\mathrm{6}}$A score RCC patients. (A–D) Violin plots of the estimated IC50 values of chemotherapy drugs including axitinib, pazopanib, sorafenib and sunitinib in high- and low-m${}^6$A score RCC patients. (E) Differential drug response analyses of six Cancer Therapeutics Response Portal (CTRP)-derived compounds in high- vs low-m${}^6$A score RCC groups. Comparisons between two groups were analyzed with Wilcoxon test, ***p 0.001. (F) Differential drug response analyses of ten PRISM-derived compounds in high- vs low-m${}^6$A score RCC groups. Comparisons between two groups were analyzed with Wilcoxon test, ***p 0.001. (G,H) Associations between m${}^6$A score and CRISPR-Cas9 essentiality screens (CERES) score and protein expression of druggable targets.