Identification of New m6A Methylation Modification Patterns and Tumor Microenvironment Infiltration Landscape that Predict Clinical Outcomes for Papillary Renal Cell Carcinoma Patients

N6-methyladenosine (m⁶A) is the product of the most prevalent mRNA modification in eukaryotic cells. Accumulating evidence shows that tumor microenvironment (TME) plays a pivotal role in tumor development. However, the underlying relationship between m⁶A modification and the TME of a papillary renal cell carcinoma (PRCC) is still unclear. To investigate the relationship between m⁶A modification and prognosis and immunotherapeutic efficacy for PRCC, we looked for distinct m⁶A modification patterns based on 23 m⁶A-related genes. Next, the correlation between m⁶A modification patterns and TME-related characteristics was investigated. Then, the intersected differentially expressed genes were selected and the scoring system, denoted as m⁶A score, was established to evaluate m⁶A modification, prognosis, and immunotherapeutic efficacy. In this study, three distinct m⁶A expression clusters were identified. Based on the results of immune cell infiltration analysis and functional analysis, carcinogenic pathways, TME-related immune cells, and pathways were identified as well. More importantly, the established m⁶A score showed good value in predicting clinical outcomes according to results using external cohorts. Specifically, PRCC patients with low m⁶A score value showed better survival, immunotherapeutic response, and higher tumor mutation burden. Furthermore, immunohistochemistry using PRCC clinical samples from our medical center was carried out and verified our results. In conclusion, this study highlights the underlying correlation between m⁶A modification and the immune landscape and, hence, enhances our understanding of the TME and improved the therapeutic outlook for PRCC patients.

Introduction

Kidney cancer is a heterogenous disease for which several subtypes with different genetic and morphologic characteristics are identified. Renal cell carcinoma (RCC) accounts for the vast majority of histological types of kidney cancer with clear cell renal cell carcinoma (ccRCC) making up 70%–80% and papillary renal cell carcinoma (PRCC) 15%–20% of RCCs (Linehan et al., 2016; Barata and Rini, 2017; Vuong et al., 2019). Although most cases of PRCC are indolent with limited risk of mortality, the overall prognosis for PRCC remains limited (Steffens et al., 19902012).

The tumor microenvironment TME is a cellular environment in which tumor cells and other nonmalignant cells exist, and it is composed of various immune cells and related materials, including lymphocytes, fibroblasts, stromal cells, blood vessels, and so on (Wu and Dai, 2017). The TME acts as the soil of tumor cells, and the great impact of TME on tumorigenesis and tumor immunotherapy has become increasingly evident (Li et al., 2021). In an abnormal TME, immune cells become significantly remodeled, which affects their normal functions, such as proliferation, migration, and differentiation (Binnewies et al., 2018). Therefore, immunosuppression is the essential characteristic of TME. Currently, RCC tumors are considered to be immunogenic, and many studies find that various immune cells could infiltrate into RCC TMEs. However, these immune cells block the effective antitumor responses. Owing to the immunosuppressed state of RCC tumors and the immune-tolerance of TMEs, the response of RCC to immune checkpoint inhibitors (ICIs) is unsatisfactory (Syn et al., 2017).

Due to the advances in RNA sequencing, N6-methyladenosine (m⁶A), the product of the most common type of mRNA modification in eukaryotic cells, has garnered great interest (Qi et al., 2016; Ke et al., 2017). The m⁶A modification is regulated by three types of molecules, known as “writer,” “eraser,” and “reader” molecules (Yang et al., 2018). It is reported that m⁶A modification plays multifaceted roles in tumor development and metastasis (Xiao et al., 2018). Various research investigation indicates that abnormal m⁶A modification occurs in most immune cells, including dendritic cells, regulatory T cells, macrophages, CD4⁺ T cells, and CD8⁺ T cells, and results in tumor escape or immune disorder (Chen et al., 2018; Han et al., 2019; Li et al., 2021). However, it is still unclear whether m⁶A modification in diverse immune cells in the TME is responsible for tumor progression and the effectiveness of ICIs. Therefore, it is essential to determine the potential effects of m⁶A modification on the TME and to explore its clinic value as a new therapeutic tool for treatment of PRCC.

Materials and Methods

Data Collection and Processing

The expression data and clinical information for kidney renal papillary cell carcinoma (KIRP) were downloaded directly from the Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/), Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/), and the Cancer Immunome Atlas (TCIA) (https://tcia.at/home). Specific data from 289 KIRP patients and 32 tumor-free patients were obtained from these databases. Copy number variation (CNV) and somatic mutation data were downloaded from TCGA as well. Samples without survival data were removed. The “limma” package was used to normalize gene expression data and transform fragments per kilobase per million (FPKM) values to transcripts per kilobase per million (TPM) value. R (R version 4.0.1) was used to extract and analyze expression data and clinical information. After conducting a comprehensive literature review (Zhang et al., 2020; Gu et al., 2021; Zhong et al., 2021), we identified 23 m⁶A regulators, including METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, RBM15B, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, RBMX, FTO, and ALKBH5, representing m⁶A writers, readers, and erasers.

Identification of Differentially Expressed Genes and Functional Analysis

The “limma” and “ggplot2” packages were used to assess and visualize the differentially expressed genes (DEGs) in KIRP samples and nontumor tissues. Difference with adjust p < .01 were considered to be significant. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed through the “clusterProfiler” package. To determine the differences in biological processes between various m⁶A expression clusters, specifically to estimate the variation in biological processes, gene set variation analysis (GSVA) was conducted by using the “GSVA” package (Hänzelmann et al., 2013). We utilized the gene set “c2.cp.kegg.v6.2-symbols” from the MSigDB database (Liberzon et al., 2011). Here, adjusted p < .05 was considered as the threshold.

Estimation of TME Immune Cell Infiltration and Tumor Mutation Burden

Single-sample gene-set enrichment analysis (ssGSEA) was used to quantify the level of immune infiltration into the PRCC TME (Barbie et al., 2009; Charoentong et al., 2017). The relevant gene set, which marks various TME-infiltrated immune cell subtypes, was collected from previous studies (Barbie et al., 2009; Charoentong et al., 2017). The ssGSEA scores represented the enrichment of different immune cell subtypes in each sample. Tumor mutation burden (TMB) was analyzed with the KIRP somatic mutation data by using the “maftools” R package (Chen and Mellman, 2017). Two TMB sets (high and low TMB) were constructed by using an optimal cutoff value of TMB. We evaluated the difference between the m6Ascore values of two TMB sets.

Unsupervised Clustering and the Construction of an m⁶A Regulators Model

Owing to relatively small sizes of the KIRP data sets in the Gene-Expression Omnibus (GEO) database, we used the GSE2748 cohort and TCGA KIRP data set to perform unsupervised clustering analysis with the “ConsensusClusterPlus” package (Wilkerson and Hayes, 2010). Here, 1000 repetitions were performed. The expression data of 23 m⁶A genes were extracted from GSE2748. The clustering analysis was performed to classify the KIRP samples into distinct m⁶A expression clusters based on the expression of 23 m⁶A regulators.

To quantify the m⁶A expression cluster of each KIRP sample, the m⁶A score was applied and established as follows. First, we identified intersected DEGs from the constructed m⁶A expression clusters. All KIRP patients were divided into diverse groups via unsupervised clustering analysis. Then, the univariate Cox regression analysis was utilized to assess the prognosis of each selected gene. p < .05 was considered as the significance criterion. After extracting the prognosis-related regulators, we applied principal component analysis to establish the m⁶A gene model, and the principal components 1 and 2 were selected as signature scores. Finally, m⁶A score was calculated using following formula: m⁶A score = Σ (PC1_i + PC2_i) (where i is the expression of the selected m⁶A related DEGs from the m⁶A expression cluster) (Sotiriou et al., 2006; Zeng et al., 2019).

Genomic and Clinical Data for ICI Therapy

Then, we investigated whether the established m⁶A expression cluster could predict the response of PRCC to ICI therapy based on two immunotherapy cohorts. After a comprehensive search for gene expression data and complete clinical information of patients treated with ICIs, we finally included two related cohorts. The first cohort involved metastatic melanoma patients treated with the anti-PD-1 drug (pembrolizumab) from the GEO database (GSE78220). Moreover, genomic and clinical data for mTOR inhibitor (everolimus) therapy was downloaded from the Supplementary File appended to published study (Barbie et al., 2009; Charoentong et al., 2017). All raw expression data were normalized using the “limma” package and transformed into the more comparable TPM value.

Immunohistochemistry

Five pairs of PRCC and adjacent normal tissues were collected from May 2021 to October 2021 from Shandong Provincial Hospital affiliated with Shandong First University. The study was approved by the Ethics Committee of Shandong Provincial Hospital (Approval No. SWYX: NO. 2021-491). IHC was performed according to published method (Wang et al., 2020). All samples were incubated with rabbit polyclonal anti-CD8 (ab101500), anti-CD69 (ab233396), anti-CD163 (ab182422), anti-YTHDF1 (ab252346), anti-YTHDF2 (ab220163), anti-YTHDF3 (ab220161), anti-ZC3H13 (IHC0104123), anti-HNRNPA2B1 (ab31645), and anti-IGFBP2 (ab188200) antibodies overnight at 4°C and then washed. Two pathologists independently assessed the IHC slides.

Statistical Analysis

The Kruskal–Wallis test was used to estimate the significance of differences between values of three or more groups. Spearman’s correlation analysis was applied to calculate the correlation coefficient between number of TME-infiltrated immune cells and the expression level of m⁶A regulators. We employed the “survminer” package to determine the optimal cutoff value. Based on the optimal cutoff point, all PRCC patients were grouped into high or low m⁶A score sets. Then, the Kaplan–Meier analysis with a log-rank test was conducted to test the prognosis of patients. The mutation landscape of KIRP cohorts was depicted by using the “maftools” package (Mayakonda et al., 2018). Statistical analysis was performed with R packages (version 4.0.1). A two-tailed p < .05 was considered to be significant.

Results

Genetic Variation and Clinical Relevance of m⁶A Genes in PRCC

Based on the transcriptomic profiles of 23 m⁶A regulators, we investigated the expression pattern of all m⁶A regulators in PRCC and normal samples from TCGA (Figure 1A). Then, we integrated CNV as well as somatic mutations and illustrated the prevalence of alteration of m⁶A genes in PRCC. Only 22 of 281 samples (7.83%) showed m⁶A regulator mutations. Specifically, 8 out of 23 m⁶A regulators experienced mutations (Figure 1B). Afterward, we investigated the CNV frequency of 23 m⁶A genes, which identified that most CNV alterations in 23 genes were focused on the CNV deletion (Figure 1C). Moreover, we determined the locations of the CNV alteration on human chromosomes as well (Figure 1D). These results indicate that genetic variation commonly occurs in PRCC cells and is heterogeneous between PRCC and normal tissues, exhibiting the potential role for the aberrant expression of m⁶A genes in tumorigenesis and development as well as progression. Finally, when investigating the potential clinical relevance of 23 m⁶A regulators, we found that three types of m⁶A regulators were positively correlated with patient prognosis and interacted with each other (Figure 1E and Supplementary Figure S1A). In addition, most of the genes were indicated to be risk factors for overall survival (OS) of PRCC patients; only YTHDC1, ALKBH5, FTO, RBM15B, METTL14, and METTL16 were out.

FIGURE 1

FIGURE 1. Characteristics of genetic and variation of m⁶A genes in pRCC. (A) Transcriptomic profile of 23 m6A genes among normal and PRCC tumor tissues. The upper and lower ends of the boxes indicate interquartile range (IQR) of values. The lines in the boxes represent each the median value. The asterisks represented the p value. (*p < .05; **p < .01; ***p < .001) (B) Mutation frequency of 23 m⁶A genes in in 289 PRCC patients. The bar plot at the upper part of the figure shows TMB. Numbers on the right part of the figure represent the mutation frequency. (C) CNV alteration frequency of 23 m⁶A genes in the 289-patient PRCC cohort. The height of each column indicates the alteration frequency. The red dot represents the amplification frequency; the green dot represents the deletion frequency. (D) Location of CNV alterations of 23 m⁶A genes on human chromosomes in PRCC cohort. (E) Interaction network of 23 m⁶A genes in the PRCC cohort. The size of each circle is indicative of the magnitude of the effect on the prognosis. The green dots mean the protective factors, and the purple dots represent the risk factors. The thickness of each line indicates the degree of correlation between each gene.

We also determined whether genetic variations of “writer,” “reader,” and “eraser” genes were associated with the expression other m⁶A regulators’ (Supplementary Figures S1B–L). The results demonstrate that only YTHDF1 was upregulated in METTL14 mutated PRCC samples while other m⁶A genes highly expressed in wild-type ALKBH5, HNRNPC, METTL14, YTHDC1, and YTHDC2.

Different m⁶A Modification Patterns Mediated by 23 m⁶A Genes and Its Clinical Relevance

Based on the expression levels of the 23 m⁶A genes, we classified the PRCC patients by carrying out unsupervised clustering analysis (Supplementary Figures S2A–E). We finally identified three patterns, termed as m⁶A expression clusters A, B and C, which included 56 cases in m⁶A expression cluster A, 128 cases in m⁶A expression cluster B, and 127 cases in m⁶A expression cluster C (Figure 2A). Then, we determined the prognostic values of the three m⁶A modification patterns. According to this analysis, m⁶A expression cluster A showed the most favorable survival (Figure 2B). After combing the TCGA and GEO data sets for comprehensive clinical data from PRCC patients, we made a heat map to visualize the correlation between the three m⁶A expression clusters and clinical characteristics. As shown in the Figure 2A, m⁶A expression cluster C was associated with poor prognosis and enriched in metastatic tumors as well as being associated with patient old age. By comparison, m⁶A expression clusters A and B showed relatively better prognose. We also noted that 23 m⁶A-related genes had relatively high expression levels in m⁶A expression cluster C, followed by m⁶A expression clusters B and A (Figure 2C).

FIGURE 2

FIGURE 2. Distinct m⁶A modification patterns and biological features of each cluster. (A) Unsupervised clustering of 23 m⁶A genes in TCGA PRCC and GSE2748 data sets. The m⁶A expression cluster, tumor M stage, survival status, gender, and age were used as annotations. Red means high expression of m⁶A genes, and blue represents low expression. The M means metastasis. (B) Kaplan–Meier survival analysis showing the OS for the three m⁶A expression clusters based on TCGA PRCC and GSE2748 data sets. A log-rank test was performed. (C) Expression pattern of 23 m⁶A genes among three m⁶A expression clusters. The upper and lower ends of the boxes mean interquartile range (IQR) of the values. The lines in the boxes represent each the median value. The asterisks represent the p value. (*p < .05; **p < .01; ***p < .001) (D,E) GSVA enrichment analysis depicting the activation states of the biological processes in three m⁶A expression clusters. In the heat map, red means activated pathways, and blue means inhibited pathways. (D) m⁶A expression clusters B vs. A (E) m⁶A expression clusters B vs. C.

Biological and TME Cell Infiltration Characteristics in Three m⁶A Modification Patterns

To investigate the biological processes associated with the three types of m⁶A modification patterns, we performed a GSVA analysis. The m⁶A expression cluster A was found to be associated with immune activation processes, such as complement and coagulation cascades. The m⁶A expression cluster B was found to be associated with oncogenic and stromal signaling pathways, including mTOR signaling pathways, ERBB signaling pathways, and adherens junction. The m⁶A expression cluster C was also found to be related with immune-related pathways, such as the Notch signaling pathway (Figures 2D,E). Then, we explored the TME cell infiltration for the different m⁶A expression clusters. The ssGSEA analysis presented that activated CD8⁺ T cells, myeloid-derived suppressor cells, and several innate immune cells, such as macrophages and monocytes, were enriched in m⁶A expression cluster A (Figure 3A). Moreover, m⁶A expression cluster C was associated with natural killer cells, plasmacytoid dendritic cells, and type 2 T helper cells. Afterward, we determined the proportion of immune cells in the three m⁶A expression clusters by using the CIBERSORT algorithm (Figure 3B). However, a significant difference between the different immune cells was not observed. Finally, we used principal component analysis (PCA), which verified significant differences between the three distinct clusters of PRCC patients (Figure 3C).

FIGURE 3

FIGURE 3. Characteristics of TME cell infiltration in diverse m⁶A expression clusters and the construction of m⁶A gene signatures. (A) Abundance of TME-related immune cells in three m⁶A patterns. The histogram shows the expression difference of 23 kinds of immune cells between m⁶A expression cluster A, B, and C. The upper and lower ends of the boxes mean interquartile range (IQR) of the values. The lines in the boxes represent each the median value. The asterisks represent the p value. (*p < .05; **p < .01; ***p < .001) (B) Proportion of different TME-related immune cells in the three m⁶A patterns as analyzed by CIBERSORT. (C) Principal component analysis of PRCC patients in three m⁶A expression clusters, which indicates a remarkable difference between the different modification patterns. (D,E) Functional analysis for m⁶A-related genes. (D) GO enrichment analysis. (E) KEGG enrichment analysis. MF means molecular function, CC means cellular component, and BP means biological process. (F) Principal component analysis of PRCC patients for the three m⁶A-based gene expression clusters. The left part of each figure indicates the biological functions and signaling pathways. The degree of enrichment is represented by the color depth of each bar plot.

Model and Biological Characteristics of the m⁶A Regulators

To further describe the features of the three m⁶A expression clusters, we identified 4780 intersected m⁶A DEGs among the three clusters (Supplementary Figure S2E). Afterward, we analyzed these phenotype-related genes by carrying out KEGG and GO enrichment analyses. The GO analysis revealed a significant enrichment (FDR <0.01) of the methyltransferase complex, RNA methyltransferase activity, and activation of innate immune response (Figure 3D and Supplementary Table S1). The KEGG pathway analysis also indicated that RCC, PD-L1 expression, and the PD-1 checkpoint pathway in cancer were enriched in these selected m⁶A DEGs (Figure 3E and Supplementary Table S2). The above analysis further confirmed the pivotal role played by m⁶A modification in immune regulation as well as RCC. Next, univariate Cox regression analysis was carried out to determine the prognosis-related m⁶A genes. Here, 1285 prognosis-related m⁶A regulators were extracted for unsupervised clustering analysis. With the optimal k = 3, three genomic clusters were constructed and named m⁶A-based gene expression clusters A–C (Supplementary Figures S3A–E). A PCA analysis found difference between these three m⁶A-based gene expression clusters as well (Figure 3F). Once again, these results confirmed that diverse m⁶A modification patterns occurred for PRCC.

To determine the clinical relevance of these clusters, we evaluated the healthy status among the three m⁶A-based gene expression clusters. The m⁶A-based gene cluster C showed a worse prognosis than did m⁶A-based gene expression clusters A and B (Figure 4A). As shown in Figure 4B, m⁶A-based gene expression cluster C was mainly enriched in metastatic tumors. However, the other clusters were related with alive status as well as nonmetastatic tumor (Figure 4B). The results of the differential analysis of the three clusters validated the pattern of m⁶Agene signatures as well (Figure 4C).

FIGURE 4

FIGURE 4. Characteristics of diverse m⁶A-based gene expression clusters. (A) Kaplan–Meier survival analysis showed the OS for the three m⁶A-based gene expression clusters based on TCGA PRCC and GSE2748 data set with PRCC. (B) Unsupervised clustering of the intersected m⁶A phenotype-related genes in PRCC, which classifies patients into several clusters, termed m⁶A-based gene expression clusters. The m⁶A cluster, tumor M stage, survival status, gender, and age are used as annotations. Red means high expression of m⁶A genes, and blue represents low expression. The M means metastasis. (C) Expression pattern of 23 m⁶A genes for the three m⁶A-based gene expression clusters. The histogram indicates the expression level of 23 m⁶A genes between m⁶A-based gene cluster A, B, and C. The upper and lower ends of the boxes mean each the interquartile range (IQR) of values. The lines in the boxes represent median value. The asterisks represented the p value. (*p < .05; **p < .01; ***p < .001). (D) Alluvial diagram displaying the differences in m⁶A expression clusters, m⁶A-based gene expression clusters, and m⁶Ascore. (E) Spearman analyses of the correlations between m⁶Ascore and biological characteristics in the PRCC cohort. (F,G) Differential analysis of m⁶A score values among (F) m⁶A expression clusters and (G) m⁶A-based gene expression clusters in TCGA PRCC and GSE2748 data sets. (H) K–M analyses for the OS of PRCC patients in high and low m⁶A score groups.

Evaluation of the m⁶A Modification Patterns Among the m⁶A Regulator Signatures

We employed m⁶A score (a scoring methodology) to quantify and evaluate m⁶A modification patterns. Alterations of each of the PRCC patient’s attributes were visualized by producing and inspecting an alluvial diagram. The results suggest that most of the PRCC samples showing the m6A-based gene expression cluster C were marked with a high m⁶A score and showed poor patient survival (Figure 4D). Then, we assessed the correlations between m⁶A score and biological processes. The m⁶A score was only positively associated with processes involving type 2 T helper cells but negatively correlated with processes involving other immune cells (Figure 4E). Significant differences in the m⁶A score were observed between the three m⁶A-based gene expression clusters as well as between the m⁶A expression clusters. Both of these results presented that m6A expression cluster C and m6A-based gene expression cluster C have the highest m⁶Ascore (Figures 4F,G). Afterward, PRCC patients were divided into two distinct groups with an optimal cutoff value. As shown in Figure 4H, patients with low m⁶A scores showed relatively good survival compared with the high m⁶A score group.

The m⁶A Modification Model in the Role of Tumor Somatic Mutation and Immunotherapy

We also analyzed and visualized the somatic mutation profiles of PRCC patients of the high and low m⁶A score groups by using the “maftools” package. Compared with the high m⁶A score set, the low m⁶A score group showed a higher percentage of somatic mutations (Figures 5A,B). A previous study shows an association of high TMB with better survival for most cancers (Xie et al., 2020). Still, a high TMB could improve the prognosis for patients treated with ICIs (Samstein et al., 2019). Considering the significant role of TMB, we tested its prognosis value for PRCC. As observed in the survival plot, the high-TMB set presented improved survival (Figure 5C). Moreover, we found the worst survival for the PRCC patients with both a low-TMB and high m⁶A score (Figure 5D). The above outcome implies that TMB as well as m⁶A score could potentially be used as predictive biomarkers.

FIGURE 5

FIGURE 5. The changes of somatic mutations among distinct m⁶A score groups. (A,B) Waterfall plot showing The changes of somatic mutations in (A) the low-m6A score set and (B) high-m6A score group. The bar plot at the upper part of the figure shows TMB. The numbers on the right part of the figure represent the mutation frequency. (C) Survival analysis for the OS of patients in high- and low-TMB groups. (D) K-M analysis for PRCC patients stratified by both m⁶A score and TMB. MSI-H, high microsatellite instability; MSI-L, low microsatellite instability.

Next, we interrogated the clinical value of the m⁶A modification model in immunotherapy (including PD-1 blockade and mTOR inhibitor). In the PD-1 blockade cohort (GSE78220), patients with low m⁶A scores showed improved overall survival (OS) (Figure 6A). In addition, in the anti-mTOR group, there was a significant difference in OS as well as progression free survival (PFS) between low and high m⁶A score groups. The therapeutic advantages of the mTOR inhibitor was observed in the low m⁶A score group (Figures 6B,C). Moreover, in light of unsatisfactory outcomes from tumor therapy, we queried whether m⁶A score could affect the therapeutic efficacy. The poor outcome of overall response rate and clinical benefit was correlated with high m⁶A score (Figures 6D,E). Finally, we used the m⁶A score to predict the reaction to immunotherapy efficacy. After downloading the immunotherapy fraction data from the Cancer Immunome Database (TCIA), we compared the predictive abilities of the m⁶A scores of the two m⁶A score groups. Patients with low m⁶A score values showed significantly better reactions to anti-CTLA-4 and anti-PD-1 therapy (Figures 6F–I).

FIGURE 6

FIGURE 6. The role of distinct m⁶A modification patterns in immunotherapy. (A–C) Results showing the associations of m6A score was negatively associated with OS and PFS following (A,B) anti-PD-1 therapy or (C) use of mTOR inhibitors. Negative associations were observed in both cases. (D,E) Proportions of PRCC patients with (D) an immunotherapy response and (E) clinical benefit in the two m6Ascore sets. (F–I) Relationship between m6A score and immunotherapeutic response under PD-1 and CTLA-4 expressions: (F) negative PD-1 and CTLA-4. (G) positive PD-1 and negative CTLA-4. (H) positive CTLA-4 and negative PD-1. (I) positive PD-1 and CTLA-4.

Biological Validation of Significant m⁶A Regulators and Immune Cell Markers

The robustness of m⁶A regulators as biomarkers was verified using primary PRCC clinical samples from the Shandong Provincial Hospital affiliated with Shandong First Medical University. We selected six m⁶A genes from the DEGs and five immune cell markers for the following validation. The IHC images acquired of immune cell markers showed weak staining for CD8, CD69, and CD163 in normal renal tissue (Figures 7G–I). Tumor tissue staining of YTHDF1 and HNRNPA2B1 showed moderate staining in the nucleus, and negative staining was observed in the normal tissues (Figures 7A,D). In normal kidney samples, moderate staining for ZC3H13 and YTHDF2 were observed in the nucleus. Regarding the YTHDF3 and IGFBP2, strong staining was positive on the cytoplasm (Figures 7B,C,E,F). However, weak staining patterns for ZC3H13, YTHDF2, YTHDF3, and IGFBP2 were observed in PRCC tissues (Figures 7B,C,E,F). These unique IHC staining patterns further confirmed the above results and illustrated that these selected m⁶A regulators could be used to predict clinical outcomes.

FIGURE 7

FIGURE 7. Representative IHC images of significant m⁶A regulators and immune cell markers in PRCC samples. (A–F) IHC patterns for selected m⁶A regulators in normal and tumor samples. (G–I) IHC patterns for three immune cell markers in normal and tumor samples. Bar, 50 and 200 μm.

Discussion

The m⁶A modification plays a pivotal role in tumorigenesis, tumor development, progression, and prognosis (He et al., 2019). Previous studies show the m⁶A modification displaying dual suppressive and promotive functions in various tumors (He et al., 2018; Wang et al., 2018). However, there are few studies of the m⁶A modification for RCC (especially PRCC) initiation, progression, and therapy. The TME is a potential regulator of cancer progression and a source of therapeutic targets. In the complex TME, immune and stromal cells play significant roles in cancer development (Quail and Joyce, 2013; Ho et al., 2020). Currently, knowledge of the kidney TME is restricted to only a few different tumor types and lacks comprehensive analysis. Therefore, in this study, we focused our attention on the role of m⁶A modification in the TME of PRCC and aimed to unravel the potential functions of this modification and contribute to obtaining a deeper understanding of antitumor immune effects of the TME in PRCC.

CNV is one of the most important somatic aberrations in cancer, and several studies find significant associations between CNVs and cancers (Speleman et al., 2008; Shlien and Malkin, 2009; Beroukhim et al., 2010). Based on 23 m⁶A genes and PRCC copy-number profiles, we explored the alteration of m⁶A genes in PRCC. The mutations of the m⁶A regulators occurred relatively infrequently in PRCC, but CNV deletion was a common event. Then, on the basis of clustering analysis, we identified three different m⁶A expression clusters in PRCC. In 2017, Chen DS et al. proposed three types of cancer-immune phenotypes, namely, immune-inflamed, immune-excluded, and immune-desert phenotypes (Speleman et al., 2008; Shlien and Malkin, 2009; Beroukhim et al., 2010). The immune-inflamed phenotype is characterized by the presence of CD4⁺ T, CD8⁺ T, myeloid, and monocytc cells in the TME, which is positioned near the tumor cells (Herbst et al., 2014; Turley et al., 2015). The immune-excluded phenotype also involves the presence of many immune cells, but with these cell, located mainly surrounding the stroma instead of the nest of the tumor (Joyce and Fearon, 2015; Hegde et al., 2016). The immune-desert phenotype presents a paucity of CD8⁺ T cells in both tumor parenchyma and stroma with this paucity being a feature of a noninflamed TME (Gajewski et al., 2013; Kim and Chen, 2016). In our current study, we found an enrichment of activated CD8⁺ T cells, myeloid-derived suppressor cells, macrophages, and monocytes in m⁶A expression cluster A, an association of the m⁶A expression cluster B with adherens junction, and m⁶A expression cluster C showing the presence of natural killer and plasmacytoid dendritic cells. Due to the presence of CD8 expressing T cells and other myeloid cells as well as monocytes, the m⁶A expression cluster A showed improved survival.

Then, we identified the intersected DEGs between diverse m⁶A expression clusters and assessed the potential biological functions of these genes and the pathways used by them. Our results show a significant enrichment of these DEGs in m⁶A-, immune- and immunotherapy-related biological functions and pathways. Moreover, we chose T cell (CD8, CD69) and macrophage markers (CD163) as well as differentially expressed m⁶A regulators to validate the clinical application using primary PRCC samples from our hospital, and the results further confirm the prognostic value in clinical application. To limit the individual heterogeneity, we utilized m⁶A score to quantify and evaluate m⁶A modification patterns. Similar to the results of previous research, the m⁶A expression cluster C and m⁶A expression cluster A in the current work presented, respectively, the highest and lowest m⁶A score in PRCC. The K-M survival curve illustrates a better OS and better prognosis associated with m⁶A-based gene expression cluster A than with m⁶A-based gene expression cluster C. These results suggest that the m⁶A scoring system could be applied to determine distinct immune phenotypes and m⁶A modification patterns.

Somatic mutation was detected between high- and low-m⁶A score groups as well. The low m⁶A score group had a high TMB with high TMB associated with better survival for PRCC patients. A similar trend was found in studies involving melanoma and osteosarcoma (Aoude et al., 2020; Xie et al., 2020). Still, a high TMB appears to indicate a better prognosis for patients receiving ICIs for treating various types of tumors (Snyder et al., 2014; Rizvi et al., 2015; Van Allen et al., 2015; Rosenberg et al., 2016). These findings suggest better immunotherapeutic outcomes for the low m⁶A score group than for the high m⁶A score group. In light of the disappointing outcomes from immunotherapy (including anti-PD-1 therapy and mTOR inhibitors) to date (Larkin et al., 2015; Postow et al., 2015; Rotte et al., 2018; de Vries-Brilland et al., 2020), we sought to determine whether m⁶A score could serve as a biomarker to stratify patients with different levels of immune-responsiveness to tumors. By utilizing GSE78220 (PD-1 blockade cohort) and the anti-mTOR group (Hugo et al., 2016; Braun et al., 2020), we showed an association between a low m⁶A score and improve OS and PFS time. Thus, distinct m⁶A modification patterns may impact the efficacy of immunotherapy, and m⁶A score has potential clinical value in evaluating the efficacy of therapeutic.

To improve the outcome for PRCC patients, access to accurate and efficient biomarkers is indispensable. Therefore, we investigated the TME and m⁶A-related genes to reveal the associated immune cells and molecular mechanism as well as clinical value. This investigation suggests that diverse m⁶A modification patterns could affect the complexity of the PRCC TME. Moreover, the established m⁶A score was indicated by our results to have great potential as a predictive indicator to assess the distinct m⁶A modification patterns and prognose of PRCC patients. More importantly, given the high variety of responses to immunotherapy, the m⁶A score may be utilized to evaluate how tumors might react to being exposed to an immunotherapy (including anti-PD-1 therapy and mTOR inhibitors). We do note that the relatively small number of PRCC patients receiving immunotherapy may affect the predictive ability of m⁶A score. Therefore, in future investigations, expression data and clinical information from our medical center will be collected. Further experiments in vivo and in vitro will also be implemented to confirm the molecular mechanism of m⁶A-related regulators in the PRCC TME. Nevertheless, the study we carried out has enhanced our understanding of TME characteristics and improved the therapeutic landscape for PRCC patients.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee of Shandong Provincial Hospital. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

BZ and WH contributed to the study conception; FC, BZ, JW, ZY, YZ, and ZC conducted the data analysis and were responsible for writing the first draft of the paper. WH, FC, and ZN revised the paper; and all authors read and approved the final version of the manuscript.

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.818194/full#supplementary-material

Abbreviations

ccRCC, clear cell renal cell carcinoma; CNV, Copy number variation; DEGs, Differentially expressed genes; GEO, Gene-Expression Omnibus; GSVA, Gene set variation analysis; ICI, Immunological checkpoint inhibitor; m⁶A, N6-methyladenosine; PCA, Principal component analysis; PRCC, papillary renal cell carcinoma; ssGSEA, Single-sample gene-set enrichment analysis; TCGA, The Cancer Genome Atlas; TCIA, The Cancer Immunome Atlas; TMB, Tumor mutation burden; TME, tumor microenvironment.

References

Aoude, L. G., Bonazzi, V. F., Brosda, S., Patel, K., Koufariotis, L. T., Oey, H., et al. (2020). Pathogenic Germline Variants Are Associated with Poor Survival in Stage III/IV Melanoma Patients. Sci. Rep. 10 (1), 17687. doi:10.1038/s41598-020-74956-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Barata, P. C., and Rini, B. I. (2017). Treatment of Renal Cell Carcinoma: Current Status and Future Directions. CA: a Cancer J. clinicians 67 (6), 507–524. doi:10.3322/caac.21411

CrossRef Full Text | Google Scholar

Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA Interference Reveals that Oncogenic KRAS-Driven Cancers Require TBK1. Nature 462 (7269), 108–112. doi:10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Beroukhim, R., Mermel, C. H., Porter, D., Wei, G., Raychaudhuri, S., Donovan, J., et al. (2010). The Landscape of Somatic Copy-Number Alteration across Human Cancers. Nature 463 (7283), 899–905. doi:10.1038/nature08822

PubMed Abstract | CrossRef Full Text | Google Scholar

Binnewies, M., Roberts, E. W., Kersten, K., Chan, V., Fearon, D. F., Merad, M., et al. (2018). Understanding the Tumor Immune Microenvironment (TIME) for Effective Therapy. Nat. Med. 24 (5), 541–550. doi:10.1038/s41591-018-0014-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Braun, D. A., Hou, Y., Bakouny, Z., Ficial, M., Sant’ Angelo, M., Forman, J., et al. (2020). Interplay of Somatic Alterations and Immune Infiltration Modulates Response to PD-1 Blockade in Advanced clear Cell Renal Cell Carcinoma. Nat. Med. 26 (6), 909–918. doi:10.1038/s41591-020-0839-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cel Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D. S., and Mellman, I. (2017). Elements of Cancer Immunity and the Cancer-Immune Set point. Nature 541 (7637), 321–330. doi:10.1038/nature21349

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Wei, L., Law, C.-T., Tsang, F. H.-C., Shen, J., Cheng, C. L.-H., et al. (2018). RNA N6-Methyladenosine Methyltransferase-like 3 Promotes Liver Cancer Progression through YTHDF2-dependent Posttranscriptional Silencing of SOCS2. Hepatology 67 (6), 2254–2270. doi:10.1002/hep.29683

PubMed Abstract | CrossRef Full Text | Google Scholar

de Vries-Brilland, M., Gross-Goupil, M., Seegers, V., Boughalem, E., Beuselinck, B., Thibault, C., et al. (2020). Are Immune Checkpoint Inhibitors a Valid Option for Papillary Renal Cell Carcinoma? A Multicentre Retrospective Study, 136. Oxford, United Kingdom: European journal of cancer, 76–83.

CrossRef Full Text | Google Scholar

Gajewski, T. F., Woo, S.-R., Zha, Y., Spaapen, R., Zheng, Y., Corrales, L., et al. (2013). Cancer Immunotherapy Strategies Based on Overcoming Barriers within the Tumor Microenvironment. Curr. Opin. Immunol. 25 (2), 268–276. doi:10.1016/j.coi.2013.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, Y., Wu, X., Zhang, J., Fang, Y., Pan, Y., Shu, Y., et al. (2021). The Evolving Landscape of N6-Methyladenosine Modification in the Tumor Microenvironment. Mol. Ther. 29 (5), 1703–1715. doi:10.1016/j.ymthe.2021.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, D., Liu, J., Chen, C., Dong, L., Liu, Y., Chang, R., et al. (2019). Anti-tumour Immunity Controlled through mRNA m6A Methylation and YTHDF1 in Dendritic Cells. Nature 566 (7743), 270–274. doi:10.1038/s41586-019-0916-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

He, L., Li, H., Wu, A., Peng, Y., Shu, G., and Yin, G. (2019). Functions of N6-Methyladenosine and its Role in Cancer. Mol. Cancer 18 (1), 176. doi:10.1186/s12943-019-1109-9

PubMed Abstract | CrossRef Full Text | Google Scholar

He, L., Li, J., Wang, X., Ying, Y., Xie, H., Yan, H., et al. (2018). The Dual Role of N6‐methyladenosine Modification of RNAs Is Involved in Human Cancers. J. Cel Mol Med 22 (10), 4630–4639. doi:10.1111/jcmm.13804

CrossRef Full Text | Google Scholar

Hegde, P. S., Karanikas, V., and Evers, S. (2016). The where, the when, and the How of Immune Monitoring for Cancer Immunotherapies in the Era of Checkpoint Inhibition. Clin. Cancer Res. 22 (8), 1865–1874. doi:10.1158/1078-0432.ccr-15-1507

PubMed Abstract | CrossRef Full Text | Google Scholar

Herbst, R. S., Soria, J.-C., Kowanetz, M., Fine, G. D., Hamid, O., Gordon, M. S., et al. (2014). Predictive Correlates of Response to the Anti-PD-L1 Antibody MPDL3280A in Cancer Patients. Nature 515 (7528), 563–567. doi:10.1038/nature14011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ho, W. J., Jaffee, E. M., and Zheng, L. (2020). The Tumour Microenvironment in Pancreatic Cancer - Clinical Challenges and Opportunities. Nat. Rev. Clin. Oncol. 17 (9), 527–540. doi:10.1038/s41571-020-0363-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hugo, W., Zaretsky, J. M., Sun, L., Song, C., Moreno, B. H., Hu-Lieskovan, S., et al. (2016). Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell 165 (1), 35–44. doi:10.1016/j.cell.2016.02.065

PubMed Abstract | CrossRef Full Text | Google Scholar

Joyce, J. A., and Fearon, D. T. (2015). T Cell Exclusion, Immune Privilege, and the Tumor Microenvironment. Science 348 (6230), 74–80. doi:10.1126/science.aaa6204

PubMed Abstract | CrossRef Full Text | Google Scholar

Ke, S., Pandya-Jones, A., Saito, Y., Fak, J. J., Vågbø, C. B., Geula, S., et al. (2017). m6A mRNA Modifications Are Deposited in Nascent Pre-mRNA and Are Not Required for Splicing but Do Specify Cytoplasmic turnoverA mRNA Modifications Are Deposited in Nascent Pre-mRNA and Are Not Required for Splicing but Do Specify Cytoplasmic Turnover. Genes Dev. 31 (10), 990–1006. doi:10.1101/gad.301036.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J. M., and Chen, D. S. (2016). Immune Escape to PD-L1/pd-1 Blockade: Seven Steps to success (Or Failure). Ann. Oncol. 27 (8), 1492–1504. doi:10.1093/annonc/mdw217

PubMed Abstract | CrossRef Full Text | Google Scholar

Larkin, J., Chiarion-Sileni, V., Gonzalez, R., Grob, J. J., Cowey, C. L., Lao, C. D., et al. (2015). Combined Nivolumab and Ipilimumab or Monotherapy in Untreated Melanoma. N. Engl. J. Med. 373 (1), 23–34. doi:10.1056/nejmoa1504030

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Zha, X., and Wang, S. (2021). The Role of N6-Methyladenosine mRNA in the Tumor Microenvironment. Biochim. Biophys. Acta (Bba) - Rev. Cancer 1875 (2), 188522. doi:10.1016/j.bbcan.2021.188522

CrossRef Full Text | Google Scholar

Liberzon, A., Subramanian, A., Pinchback, R., Thorvaldsdóttir, H., Tamayo, P., and Mesirov, J. P. (2011). Molecular Signatures Database (MSigDB) 3.0. Bioinformatics 27 (12), 1739–1740. doi:10.1093/bioinformatics/btr260

PubMed Abstract | CrossRef Full Text | Google Scholar

Linehan, W. M., Linehan, W. M., Spellman, P. T., Ricketts, C. J., Creighton, C. J., Fei, S. S., et al. (2016). Comprehensive Molecular Characterization of Papillary Renal-Cell Carcinoma. N. Engl. J. Med. 374 (2), 135–145. doi:10.1056/NEJMoa1505917

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayakonda, A., Lin, D.-C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Postow, M. A., Chesney, J., Pavlick, A. C., Robert, C., Grossmann, K., McDermott, D., et al. (2015). Nivolumab and Ipilimumab versus Ipilimumab in Untreated Melanoma. N. Engl. J. Med. 372 (21), 2006–2017. doi:10.1056/nejmoa1414428

CrossRef Full Text | Google Scholar

Qi, S. T., Ma, J. Y., Wang, Z. B., Guo, L., Hou, Y., Sun, Q. Y., et al. (2016). N6-Methyladenosine Sequencing Highlights the Involvement of mRNA Methylation in Oocyte Meiotic Maturation and Embryo Development by Regulating Translation in Xenopus laevis. J. Biol. Chem. 291 (44), 23020–23026. doi:10.1074/jbc.m116.748889

PubMed Abstract | CrossRef Full Text | Google Scholar

Quail, D. F., and Joyce, J. A. (2013). Microenvironmental Regulation of Tumor Progression and Metastasis. Nat. Med. 19 (11), 1423–1437. doi:10.1038/nm.3394

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizvi, N. A., Hellmann, M. D., Snyder, A., Kvistborg, P., Makarov, V., Havel, J. J., et al. (2015). Mutational Landscape Determines Sensitivity to PD-1 Blockade in Non-small Cell Lung Cancer. Science 348 (6230), 124–128. doi:10.1126/science.aaa1348

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenberg, J. E., Hoffman-Censits, J., Powles, T., van der Heijden, M. S., Balar, A. V., Necchi, A., et al. (2016). Atezolizumab in Patients with Locally Advanced and Metastatic Urothelial Carcinoma Who Have Progressed Following Treatment with Platinum-Based Chemotherapy: a Single-Arm, Multicentre, Phase 2 Trial. The Lancet 387 (10031), 1909–1920. doi:10.1016/s0140-6736(16)00561-4

CrossRef Full Text | Google Scholar

Rotte, A., Jin, J. Y., and Lemaire, V. (2018). Mechanistic Overview of Immune Checkpoints to Support the Rational Design of Their Combinations in Cancer Immunotherapy. Ann. Oncol. 29 (1), 71–83. doi:10.1093/annonc/mdx686

PubMed Abstract | CrossRef Full Text | Google Scholar

Samstein, R. M., Lee, C.-H., Shoushtari, A. N., Hellmann, M. D., Shen, R., Janjigian, Y. Y., et al. (2019). Tumor Mutational Load Predicts Survival after Immunotherapy across Multiple Cancer Types. Nat. Genet. 51 (2), 202–206. doi:10.1038/s41588-018-0312-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Shlien, A., and Malkin, D. (2009). Copy Number Variations and Cancer. Genome Med. 1 (6), 62. doi:10.1186/gm62

PubMed Abstract | CrossRef Full Text | Google Scholar

Snyder, A., Makarov, V., Merghoub, T., Yuan, J., Zaretsky, J. M., Desrichard, A., et al. (2014). Genetic Basis for Clinical Response to CTLA-4 Blockade in Melanoma. N. Engl. J. Med. 371 (23), 2189–2199. doi:10.1056/nejmoa1406498

CrossRef Full Text | Google Scholar

Sotiriou, C., Wirapati, P., Loi, S., Harris, A., Fox, S., Smeds, J., et al. (2006). Gene Expression Profiling in Breast Cancer: Understanding the Molecular Basis of Histologic Grade to Improve Prognosis. J. Natl. Cancer Inst. 98 (4), 262–272. doi:10.1093/jnci/djj052

CrossRef Full Text | Google Scholar

Speleman, F., Kumps, C., Buysse, K., Poppe, B., Menten, B., and De Preter, K. (2008). Copy Number Alterations and Copy Number Variation in Cancer: Close Encounters of the Bad Kind. Cytogenet. Genome Res. 123 (1-4), 176–182. doi:10.1159/000184706

PubMed Abstract | CrossRef Full Text | Google Scholar

Steffens, S., Janssen, M., Roos, F. C., Becker, F., Schumacher, S., Seidel, C., et al. (19902012)., 48. Oxford, England, 2347–2352. doi:10.1016/j.ejca.2012.05.002Incidence and Long-Term Prognosis of Papillary Compared to clear Cell Renal Cell Carcinoma-Aa Multicentre StudyEur. J. Cancer15

CrossRef Full Text | Google Scholar

Syn, N. L., Teng, M. W. L., Mok, T. S. K., and Soo, R. A. (2017). De-novo and Acquired Resistance to Immune Checkpoint Targeting. Lancet Oncol. 18 (12), e731–e741. doi:10.1016/s1470-2045(17)30607-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Turley, S. J., Cremasco, V., and Astarita, J. L. (2015). Immunological Hallmarks of Stromal Cells in the Tumour Microenvironment. Nat. Rev. Immunol. 15 (11), 669–682. doi:10.1038/nri3902

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Allen, E. M., Miao, D., Schilling, B., Shukla, S. A., Blank, C., Zimmer, L., et al. (2015). Genomic Correlates of Response to CTLA-4 Blockade in Metastatic Melanoma. Science 350 (6257), 207–211. doi:10.1126/science.aad0095

PubMed Abstract | CrossRef Full Text | Google Scholar

Vuong, L., Kotecha, R. R., Voss, M. H., and Hakimi, A. A. (2019). Tumor Microenvironment Dynamics in Clear-Cell Renal Cell Carcinoma. Cancer Discov. 9 (10), 1349–1357. doi:10.1158/2159-8290.cd-19-0499

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Chen, C., Ding, Q., Zhao, Y., Wang, Z., Chen, J., et al. (2020). METTL3-mediated m6A Modification of HDGF mRNA Promotes Gastric Cancer Progression and Has Prognostic Significance. Gut 69 (7), 1193–1205. doi:10.1136/gutjnl-2019-319639

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Chai, P., Jia, R., and Jia, R. (2018). Novel Insights on m6A RNA Methylation in Tumorigenesis: a Double-Edged Sword. Mol. Cancer 17 (1), 101. doi:10.1186/s12943-018-0847-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a Class Discovery Tool with Confidence Assessments and Item Tracking. Bioinformatics (Oxford, England) 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, T., and Dai, Y. (2017). Tumor Microenvironment and Therapeutic Response. Cancer Lett. 387, 61–68. doi:10.1016/j.canlet.2016.01.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, C.-L., Zhu, S., He, M., Chen, D., Zhang, Q., Chen, Y., et al. (2018). N6-Methyladenine DNA Modification in the Human Genome. Mol. Cel. 71 (2), 306–318. doi:10.1016/j.molcel.2018.06.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, L., Yang, Y., Guo, W., Che, D., Xu, J., Sun, X., et al. (2020). The Clinical Implications of Tumor Mutational Burden in Osteosarcoma. Front. Oncol. 10, 595527. doi:10.3389/fonc.2020.595527

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Hsu, P. J., Chen, Y.-S., and Yang, Y.-G. (2018). Dynamic Transcriptomic m6A Decoration: Writers, Erasers, Readers and Functions in RNA Metabolism. Cell Res 28 (6), 616–624. doi:10.1038/s41422-018-0040-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, D., Li, M., Zhou, R., Zhang, J., Sun, H., Shi, M., et al. (2019). Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol. Res. 7 (5), 737–750. doi:10.1158/2326-6066.cir-18-0436

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). m6A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric cancerA Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, J., Liu, Z., Cai, C., Duan, X., Deng, T., and Zeng, G., (2021). m6A Modification Patterns and Tumor Immune Landscape in clear Cell Renal Carcinoma. J. Immunother. Cancer 9 (2), 646. doi:10.1136/jitc-2020-001646

CrossRef Full Text | Google Scholar