Comprehensive Evaluation of the m6A Regulator Prognostic Risk Score in the Prediction of Immunotherapy Response in Clear Cell Renal Cell Carcinoma.

Background: Clear cell renal cell carcinoma (ccRCC) is known for its high drug resistance. The tumor-immune crosstalk mediated by the epigenetic regulation of N6-methyladenosine (m6A) modification has been demonstrated in recent studies. Therefore, m6A modification-mediated immune cell infiltration characteristics may be helpful to guide immunotherapy for ccRCC.
Methods: This study comprehensively analyzed m6A modifications using the clinical parameters, single-cell RNA sequencing data, and bulk RNA sequencing data from the TCGA-ccRC cohort and 13 external validation cohorts. A series of bioinformatic approaches were applied to construct an m6A regulator prognostic risk score (MRPRS) to predict survival and immunotherapy response in ccRCC patients. Immunological characteristics, enriched pathways, and mutation were evaluated in high- and low-MRPRS groups.
Results: The expressional alteration landscape of m6A regulators was profiled in ccRCC cell clusters and tissue. The 8 regulator genes with minimal lambda were integrated to build an MRPRS, and it was positively correlated with immunotherapeutic response in extent validation cohorts. The clinicopathological features and immune infiltration characteristics could be distinguished by the high- and low-MRPRS. Moreover, the MRPRS-mediated mutation pattern has an enhanced response to immune checkpoint blockade in the ccRCC and pan-cancer cohorts. Conclusions: The proposed MRPRS is a promising biomarker to predict clinical outcomes and therapeutic responses in ccRCC patients.

Introduction

Clear cell renal cell carcinoma (ccRCC) is the most common type of renal cancer and accounts for nearly 3% of adult malignant tumors (1). Approximately 30% of patients already have advanced ccRCC or metastases when they are first diagnosed, and have missed the opportunity for surgical intervention (2). Although targeted therapy and immunotherapy have become the main adjuvant therapy for advanced ccRCC, the complete response rate and partial response rate remained low (3, 4). So far, the biomarker-based therapeutic strategies for advanced ccRCC have been missing.

N6-methyladenosine (m6A) modification is an important factor for messenger RNA (mRNA) stability, splicing, and translation (5–7). Serval m6A-sequencing studies have revealed that abnormal m6A regulatory enzymes are involved in mutagenesis, proliferation, and tumorigenesis through the dysregulation of the m6A pathway (8, 9). Recently, m6A modifications have been shown to play a role in the regulation of immune cells, such as the following: METTL3-mediated m6A modification increased the translation of certain immune transcripts and physiologically promoted the activation of dendritic cells (DCs) and DC-based T-cell responses (10), and ALKBH5 regulated m6A modification in the 3’UTR region of PD-L1 mRNA and inhibited the expansion and cytotoxicity of T cells by sustaining tumor cell PD-L1 expression (11). The potential relationship between RNA m6A dysregulation and tumor-infiltrating immune cells (TIICs) has motivated us to investigate and find the potential biomarkers for predicting immune checkpoint therapy outcomes. Herein, we systematically evaluated the m6A regulator-based risk score and its associated gene mutation with the TIICs and revealed a new predictive method that could be used to predict the immunotherapy response in ccRCC and pan-cancer.

Materials and Methods

Data Collection and Processing

The RNA sequencing (RNA-seq) transcriptome data of patients with ccRCC and the corresponding clinical data and mutation profiles were downloaded from The Cancer Genome Atlas (TCGA) database. The validation datasets (GSE53757, GSE40435, GSE29609, and E-MTAB-3267) were included for analysis from the Gene Expression Omnibus (GEO) database and the European Molecular Biology Laboratory (12–15). The relative transcriptomic and clinical data of three immunotherapeutic cohorts of patients with ccRCC were obtained from the online supplementary data (16–18). RNA-seq data of the immunotherapy cohort of bladder cancer (19) and melanoma (PRJEB23709 and phs000452) were collected for testing (20, 21). The annotated response and mutational data of patients from a discovery cohort receiving ICB treatment from 4 studies were collected and consolidated to study the relationship between mutated genes and immunotherapy (17, 22–24). The single-cell dataset of ccRCC ICB treatment was obtained from PMID33861994 (25). The information for all collected data is presented in .

Single-Cell RNA Sequencing Analysis

The association with m6A regulators was established by analyzing the genes related to the immune response in the scRNA-seq results of ccRCC (25). The CellRanger software (version 5.0.0) and STAR were used for preprocessing. Principal component analysis (PCA) was run using the “RunPCA” function on the variable genes identified, and the k-nearest neighbor graph was constructed by the “FindNeighbors” function. Uniform manifold approximation and projection (UMAP) was used to visualize single-cell transcriptional profiles and clusters. Marker genes were visualized on UMAP plots using log-normalized counts.

Cell–Cell Communication Analysis

CellPhoneDB applies an algorithm that considers only receptors and ligands with broad expression among the tested cell types, followed by calculating the likelihood of cell-type specificity of a given receptor–ligand complex with a sufficient number of permutations (26).

Selection of m6A RNA Methylation Regulators

Based on previous studies (5, 27–30), 23 m6A RNA methylation regulators, namely, ALKBH5, CBLL1, FMR1, IGF2BP1/2/3, FTO, YTHDC1/2, YTHDF1/2/3, HNRNPC, LRPPRC, METTL3/14/16, WTAP, KIAA1429, RBM15/15B, ZC3H13, and HNRNPA2B1, were used for our analysis. Immunohistochemistry (IHC) images of m6A regulators have been used in the tissue atlas and pathology atlas panels in the Human Protein Atlas. The protein and gene expression of m6A regulators in normal individuals and ccRCC patients were analyzed on University of Alabama Cancer Database (UALCAN) portal.

Immune Infiltration Analysis in RCC

CIBERSORT and MCP counter were used to transform the RNA-seq data into the proportion of TIICs. The MCP counter R package was used to evaluate the expression of nine TIICs types. CIBERSORT (https://cibersort.stanford.edu/) was used to quantify the 22 infiltrated immune cells according to normalized gene expression profiles, which included different types of B cells, T cells, NK cells, DC cells, and mast cells. As a verification method, the single-sample gene-set enrichment analysis (ssGSEA) and xCell algorithm were applied.

Construction and Validation of the m6A Gene Signature

The significant m6A RNA methylation regulators were established by the least absolute shrinkage and selection operator (LASSO) Cox regression (with the penalty parameter estimated by 20-fold cross-validation). Those regulator genes with minimal lambda were integrated to build an MRPRS, and it was developed according to the expression level using univariate Cox. The “glmnet” package was used to perform the LASSO Cox regression model analysis.

The limma R package’s empirical Bayesian approach was applied to determine differentially expressed genes (DEGs) between high and low m6A scores. The significance criteria for determining DEGs were set as the adjusted p < 0.05 and |logFC| > 1. Finally, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the ClusterProfiler R package based on these DEGs. A protein–protein interaction (PPI) network was constructed by STRING (https://string-db.org/) and evaluated using the Cytoscope software (31).

Statistical Analysis

Statistical tests were carried out using R version 4.0.4, SPSS 25.0 (IBM, NY, USA) and GraphPad Prism 8.0. The expression levels of the m6A RNA regulators were compared with the Mann–Whitney U test in ccRCC versus normal tissues. Survival curves were generated using the Kaplan–Meier method, and the difference was compared with the log-rank test. Pearson correlation coefficient was used to compare the correlation between MRPRS and gene expression values. The “oncoplot” function of the R package “maftools” was used to determine the mutation landscape of the TCGA ccRCC cohort and immunotherapeutic cohort. The high- and low-group was divided based on the optimal cut-off value calculated by the function "surv_cutpoint" in the R package "survminer". All the R package used in this study is listed in . p < 0.05 indicated statistical significance.

Results

The Expressional Alteration Landscape of m6A Regulators in ccRCC Tissues and Cell Clusters

On reviewing the literature (5, 27–30), 23 genes were found that mainly regulate m6A modification including 9 writers (METTL3, METTL14, METTL16, RBM15, RBM15B, CBLL1, ZC3H13, KIAA1429, and WTAP), 12 readers (FMR1, HNRNPA2B1, HNRNPC, YTHDF1/2/3, YTHDC1/2, IGF2BP1/2/3, and LRPPRC), and 2 erasers (FTO and ALKBH5). We utilized the bulk TCGA-ccRCC data (529 cases of ccRCC and 74 cases of normal tissues) to analyze the expression of these m6A regulators, revealing that 15 out of 23 m6A regulators were differentially expressed (; ), and this phenomenon was also found in two other GEO datasets (GSE53757 and GSE40435) (). In addition, the protein levels of these regulators were also evaluated from the IHC results (; ).

Figure 1

(A) Expression of 23 m6A RNA methylation regulators between renal cancer and normal tissues in the TCGA-ccRCC cohort. (B, C) The UMAP plot and overview of epithelial cells by the origin and cell type of the cells. (D) Composition of various epithelial cells in different immunotherapeutic responses. (E) The heatmap of marker gene expression in 7 identified epithelial cell subsets. (F) Dot plot analysis of KEGG pathway enrichment of 7 epithelial cell subsets. (G) Violin plots showing the partial expression of m6A regulators for each epithelial cell type ("*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; ns, not significant).

We next used the scRNA-seq data (25) (PMID33861994) to evaluate the expression of m6A regulators in different subsets of ccRCC cells. A total of 65,535 cells were divided into 7 cell clusters, a total of 6,539 epithelial cells from ccRCC multiple regions (Near, Far, Center, and Lymph node) and normal tissues of 6 ICB-treated and untreated patients were extracted, and 6 ccRCC cell clusters were identified based on AQP3, GPX3, CCNI, STMN1, VCAM1, and VCAN expression (; ; ). The genes with the most significant differential expression in each cell cluster were described in the heatmap (; ), and significant functional heterogeneity was found among the 7 cell clusters (). As shown in , m6A regulators also presented the expression heterogeneity between normal cell clusters and tumor cell clusters, as well as between the 6 ccRCC cell clusters. Moreover, the expression of WTAP, YTHDC1, YTHDC2, HNRNPC, and HNRNPA2B1 was significantly different between the ICB-resistant-related GPX3+ epithelial cells and ICB-response-related VCAM1+ epithelial cells (), indicating that the differential expression level of m6A regulators in ccRCC tissue and cell clusters may be related to the efficacy of immunotherapy.

Construction, Validation, and Immunotherapy Response Evaluation of the m6A Regulator Prognostic Risk Score

To systematically evaluate these differences in m6A regulators, the MRPRS was established by the LASSO Cox algorithm, the 8 regulator genes with minimal lambda were integrated to build an MRPRS () (), and it was developed according to the expression level using univariate Cox (). The ccRCC patients in the TCGA database were divided into the high-MRPRS group (N = 134) and the low-MRPRS group (N = 395) based on the optimal cut-off value calculated by the function "surv_cutpoint" in the R package "survminer", and the ccRCC patients in the high-MRPRS group had a significantly shorter overall survival time than that in the low-MRPRS group (p < 0.001; ). The prognostic value of MRPRS was also validated in an independent cohort (GSE29609, p = 0.037; ). We continued to extend the MRPRS signature to 16 other tumor types, such as esophageal carcinoma (ESCA), lung adenocarcinoma (LUAD), and pancreatic adenocarcinoma (PAAD) (). These results present that MRPRS is negatively associated with survival outcomes.

Figure 2

(A) LASSO coefficient profiles of the 23 m6A RNA methylation regulators in the TCGA-ccRCC cohort. (B) The prognostic analyses for 23 m6A RNA methylation regulators in the TCGA-ccRCC cohort using the univariate Cox regression model. (C) Kaplan–Meier analysis of patients between high- and low-MRPRS groups in the TCGA-ccRCC cohort. (D) Validation cohort of MRPRS from GSE29609. (E) Box plot of different MRPRSs in 3 epithelial cell subsets. (F) The MRPRS between response and non-response groups in PMID32472114. (G–I) Kaplan–Meier analysis of three validation cohorts of immunotherapy in ccRCC (PMID32472114, PMID29301960, and PMID32895571). (J) Heatmap of chemotherapy and targeted drug-related genes between high- and low-MRPRS groups. (K) The correlation of MRPRS and genes associated with angiogenesis. (L) The MRPRS between response and non-response groups in ccRCC with sunitinib (E-MTAB-3267). (M) Kaplan–Meier analysis between high- and low-MRPRS groups in ccRCC with sunitinib (E-MTAB-3267) (*p < 0.05; **p < 0.01; ***p < 0.001).

Next, we investigated the correlation between MRPRS and immunotherapy response in three independent ccRCC cohorts (PMID29301960, PMID32472114, and PMID32895571) (16, 18, 22), and found that the MRPRS was significantly higher in the response group than in the non-response group; the high-MRPRS group presented a markedly prolonged survival (). Moreover, the increased MRPRS in the VCAM1+ cell cluster was positively correlated with the patients who experienced complete and mixed responses (). Similar results were also obtained in the extended dataset of bladder cancer (IMvigor210) and melanoma (PRJEB23709) ().

In addition, we analyzed the expression of targeted therapy- and chemotherapy-related genes between high- and low-MRPRS groups (). Interestingly, VEGF and mTOR pathway-related genes were found to be highly expressed in the low-MRPRS group, and MRPRS was negatively correlated with the expression level of angiogenesis-related genes including PECAM, FLT1/4, VWF, and CDH5 (). In the E-MTAB-3267 cohort of ccRCC patients treated with sunitinib, the MRPRS was significantly lower in the response group than in the non-response group, and the low-MRPRS group showed a markedly prolonged survival (). Collectively, our data suggest that the patients with high MRPRS may benefit from immunotherapy and those with low MRPRS may benefit from targeted therapy.

The Clinicopathological Features and Immune Infiltration Characteristics in Distinct MRPRS

We examined the correlation between the MRPRS and the clinical parameters. No significant association was found between the MRPRS and gender and age, but significant associations in terms of TNM stages, grade, and survival status were observed (), and MRPRS is positively correlated with TNM stages and grades (p < 0.001; ). Moreover, each of the four different T stages including stage I had significantly higher MRPRS when compared with the control subjects, and MRPRS in the metastasis group was significantly higher than that in the non-metastasis group (p < 0.001, respectively; ).

Figure 3

(A) Box plot of the relationship between stage T, N, M, grade, and MRPRS. (B) The heatmap of markers on multiple immune infiltrates. (C) The MCP_counter algorithm was used to estimate the abundance of various types of immune cells between high- and low-MRPRS groups. (D) The abundance and diversity of TCR clone in high- and low-MRPRS groups. (E) Crosstalk between immune cells and epithelial cells (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).

To investigate the effects of MRPRS on the immune infiltration characteristic of ccRCC, we evaluated the expression of immunomodulators and the infiltration levels between high- and low-MRPRS groups in ccRCC, as shown in ; 5 immunomodulators (chemokine, receptor, immunostimulator, inhibitory immune checkpoint, and MHC) and the infiltration levels of 4 types of TIICs (CD8+ T cells, DC, macrophages, and Th1 cells) were positively correlated with the high-MRPRS group (p < 0.05). The MCP counter, xCell, CIBERSORT, and ssGSEA algorithm were used to calculate an immune score and to estimate the abundance of various types of immune cells. We found significantly higher estimates of Tregs, CD8+T cells, NK cells, and B cells in ccRCC with high MRPRS ( and ; ). Moreover, the VCAM1+ cell cluster presented upregulation of HLA-A and HLA-B (), and it was also the cell cluster that communicates most frequently with immune cells (; ).

We performed volcano plots based on the DEGs from the high- and low-MRPRS groups. The results of the volcano plots showed that 1,780 genes were significantly upregulated in the comparison of the high- and low-MRPRS groups (). In the PPI network from the STRING database with the Cystoscope software, we constructed a co-expression network consisting of 45 nodes and 169 edges (). These included immune-related genes, CD19 and CD79A, and membrane proteins on the surface of B cells, which participate in the proliferation and differentiation of B cells. FOXP3 and IL2RA (CD25) are the characteristic markers of Treg cells. We also found that the expression of many costimulatory factors, such as TNFSF14, TNFRSF18, and a large amount of interleukins such as IL2 and IL6, promotes T-cell proliferation and T-cell-mediated killing (). GO enrichment analysis and KEGG analysis of these signature genes revealed that these DEGs were enriched in several biological processes and pathways related to immune regulation (; ). Moreover, the number and diversity of T-cell receptors (TCRs) were higher in the high-MRPRS group than in the low-MRPRS group (p < 0.01) (). These findings suggest that the regulation of TCR gene expression may be influenced by the specific tumor cell cluster with abnormal m6A modifications.

The Landscape of Genetic Variation of MRPRS Groups in ccRCC

The somatic mutation profile between the high- and low-MRPRS groups in the TCGA-ccRCC cohort used the maftools package, and the top 10 most frequently mutated genes in each group are shown in . Notably, SETD2, TRIOBP, RYR2, ZFPM2, and ABCC6 occupy the top 5 positions among differently mutated genes between the high- and low-MRPRS group (), and a lollipop plot showed the different mutation spots of these mutated genes between two groups (). Interestingly, the mutation rate of SETD2 was 23.18% in the high-MRPRS group and 3.87% in the low-MRPRS group, and the remaining four genes were mutated only in the high-MRPRS group. In addition, the distribution of variants according to variant classification, variant type, and single-nucleotide variant (SNV) class was displayed as a cohort summary plot, and among all the genomic alterations, missense mutations were the predominant type, with C>T and C>G representing the most common SNV classes (). Somatic mutation gene interaction networks showed a high correlation between VHL and PBRM1, PBRM1 and SETD2, and TTN and MUC16 in the high MRPRS score group ().

Figure 4

(A) Waterfall plot of the distribution of mutations found in the high- and low-MRPRS groups of the TCGA-ccRCC cohort. (B) The top 5 genes of high- vs. low-MRPRS group mutation status. (C) Lollipop plot of somatic mutations in SETD2, TRIOBP, RYR2, ZFPM2, and ABCC6. (D, E) MCP_counter and ssGSEA algorithm were used to estimate the abundance of various types of immune cells in high- and low-MRPRS groups (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; ns, not significant).

We also applied the MCP counter and ssGSEA algorithm to estimate the tumor-infiltrating immune cells between the group defined by patients with at least one mutation in these five genes or without mutation. As shown in , the infiltrating immune cells of T cells, DC cells, and B cells in the mutation status group were higher than those in the non-mutation status group (p < 0.01).

The Role of the MRPRS-Mediated Mutation Pattern in Predicting the Response to Immunotherapy

We next investigated whether the MRPRS-mediated mutation pattern could predict patients’ response to immunotherapy. We constructed a pan-cancer cohort with anti-PD-1/PDL1 immunotherapy consisting of 1,959 cases based on four cohorts (17, 22–24) (), and patients with mutation exhibited a significantly clinical response to immunotherapy and markedly prolonged survival in ccRCC (). Immunotherapy represented by PD-L1 and PD-1 blockade is a breakthrough in tumor therapy. We continued to extend the potential role of MRPRS-mediated mutation pattern in predicting responses to immunotherapy in pan-cancer () and revealed that the OS and PFS in patients with mutations were significantly higher than in those without mutations (). However, the MRPRS-mediated mutation pattern had no significance in OS of either TCGA-ccRCC or TCGA-pan-cancer (); by contrast, the PFS of the mutation group was worse than that of the non-mutation group in the TCGA-ccRCC (p = 0.049, ), and the PFS in the TCGA-pan-cancer was not significant ().

Figure 5

(A) Different distribution ratio of response and non-response in the immunotherapeutic cohort of ccRCC. (B) Kaplan–Meier analysis of patients in the mutated and non-mutated groups in the immunotherapeutic cohort of ccRCC (Van_2018, Morris_2019, PMID29337640 and PMID29301960). (C) The composition of major cancer types in the immunotherapeutic cohort contains mutations of pan-cancer. (D, E) Kaplan–Meier analysis (OS and PFS) of patients in the mutation and non-mutation groups in the immunotherapeutic cohort of pan-cancer. (F) Kaplan–Meier analysis of patients in the low TMB of mutated, low TMB of non-mutated, and high TMB groups.

In addition, tumor mutation burden (TMB) may serve as a biomarker for predicting the response to ICB treatment. We next divided the pan-cancer cohort patients into three groups according to TMB and MRPRS-mediated mutations and found that the OS of patients with low TMB and mutations was significantly better than that of the patients with high TMB and the patients with low TMB and non-mutations (P < 0.0001, ).

Discussion

The TME of ccRCC is known to be highly immunosuppressive (32). In the TME, T cells are continuously exposed to antigens, which leads to the impairment of T-cell function and ultimately to a dysfunctional state called “exhaustion” (33). The use of monoclonal antibodies or small molecules to reverse T-cell exhaustion is the basic strategy of immunotherapy (34). Since the results of the Checkmate-025 study, the immunotherapy of ccRCC has been the focus of attention, and now, combined targeted and immunotherapy has become a key component of the adjuvant treatment of advanced ccRCC (35). However, the complete or mixed response rate of immunotherapy in ccRCC is still low. Relying on biomarkers to screen patients who benefit from immunotherapy and to avoid overtreatment has long been expected in clinical practice.

Increasing evidence has demonstrated that m6A modification plays an indispensable role in immunity, inflammation, and therapy resistance through various m6A regulators (36). In this study, we systematically evaluated the expression level of m6A regulators in the ccRCC tissue and cell clusters and focused on the detailed relationship between m6A modification and TME to enhance our understanding of the ccRCC-immune crosstalk. We constructed an MRPRS comprising 8 m6A regulators by the LASSO algorithm to provide reliable biomarkers able to predict the prognosis and immunotherapy efficacy. For the first time, we analyzed the MRPRS levels in ccRCC cell clusters and found that the increased MRPRS in the VCAM1+ cell cluster was positively correlated with patients who experienced complete and mixed responses. This is consistent with our finding that the positive correlation between MRPRS and immunotherapy benefits the bulk tissue datasets. It is interesting to note that the spatial localization of this immunotherapy-related ccRCC cell cluster is worthy of further investigation.

We further explored the detailed role of m6A modification in modifying immune characteristics in ccRCC. The results of the GO and KEGG pathway analyses revealed a significant enrichment of genes in immune-related pathways. GO enrichment analysis showed that these DEGs were enriched in the humoral immune response, immunoglobulin complex, and antigen binding. The results of the KEGG analysis indicated these enriched pathways such as neuroactive ligand–receptor interaction, cytokine–cytokine receptor interaction, and the calcium signaling pathway. These results indicated that DEGs in the ccRCC are enriched in immune-related genes distinguished by the MRPRS. Among these DEGs, numerous immune-related genes were found, such as CD19, CD79A, FOXP3, CXCL13, IL2, and TNFRSF13B. FOXP3 is a hallmark of regulatory T cells, CXCL13 is related to CD8 T cells, and CD79A, CD19, and TNFRSF13B are markers of neoplastic B cells. This was in accordance with results from the single-cell sequencing analysis of ccRCC (25). Moreover, these immune cells comprise the main part of tertiary lymphoid structures (TLSs), which have recently been associated with effective antitumor immune responses in cancer patients (37, 38). These findings suggest that m6A modification may influence the formation of tertiary lymphatic structures.

The patient with a high MRPRS has a poor prognosis, and this could be due to the observation that several critical inhibitor immune checkpoints were significantly highly expressed in the high-MRPRS group, which may limit cytotoxic immune cell activities in the TME, such as CD8 T cells, causing cytotoxic cells to be in an exhausted functional state (39). Therefore, patients with high MRPRS may be more sensitive to immunotherapy. Several studies have also demonstrated that inflammatory tumor phenotypes are more sensitive to ICB (40, 41). We next compared the prognostic value of the MRPRS based on ccRCC immunotherapeutic cohorts, and the high-MRPRS group presented a prolonged survival. These findings suggest that MRPRS could be used as a new predictive biomarker for immunotherapy response in ccRCC.

Furthermore, we identified 5 genes (SETD2, TRIOBP, RYR2, ZFPM2, and ABCC6) that show the most significant differences in the comparison of mutated genes between two MRPRS groups. We found that the patients with mutated genes had worse PFS outcomes than the non-mutated group, and this was consistent with the high-MRPRS group showing worse survival than those with lower MRPRS. The potential association of TMB with sensitivity to ICB is based on the hypothesis that in tumors with high TMB, there is an increased production of surface neoantigens, thus stimulating the anti-tumor immune system response (42). The TMB has been investigated in several tumor settings, mainly in NSCLC and melanoma, as a stratification marker to predict the response to immune agents, showing promising yet inconclusive results (43, 44). In contrast, it has also been reported that high TMB fails to predict immune checkpoint blockade response across all cancer types (45). Herein, we applied the prediction of immunotherapeutic efficacy with the MRPRS-mediated mutation pattern and TMB in pan-cancer cohort and found that the MRPRS-mediated mutation pattern was a better predictor of immunotherapy outcome than the TMB. The regulatory relationship between the m6A modification and the gene mutation still needs to be studied.

Consequently, we provided a new perspective on the immune characteristics and immunotherapy strategies of ccRCC. However, several limitations should be recognized. Although we analyzed immune cell characteristics in a scRNA-seq dataset, the tumor-infiltrating immune cells were obtained based on algorithms, and thus, further experimental validation in vitro/in vivo is needed. Our study was also limited by the lack of clinical datasets to verify the relationship between the MRPRS and patients receiving targeted treatment or ICB combined targeted treatment. The combination of an MRPRS-based panel with prospective clinical trials is worth carrying out in the future.

Conclusion

This study revealed a significant association between MRPRS and TIICs of ccRCC. The proposed MRPRS is a promising biomarker to predict clinical outcomes and therapeutic responses in ccRCC patients.

Data Availability Statement

All data used in this work can be acquired from the Gene-Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under the accession numbers GSE53757, GSE40435 and GSE29609, The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/), The European Molecular Biology Laboratory (https://www.embl.org/), UALCAN (http://ualcan.path.uab.edu/), HumanProtein Atlas (HPA, https://www.proteinatlas.org/) and https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3267/. The authors would like to thank the above database for the data provided.

Author Contributions

AZ and ZZ designed this work. MY, XL, HX, SS, FW, DC, and ZW integrated and analyzed the data. AZ, MY, and XL wrote this manuscript. AZ, ZZ, ZW, and GL edited and revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (Reference number: 81402117) and the Qianjiang Talent Project of Zhejiang Province (Reference number: QJD1602025)

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.