Application of m6A and TME in Predicting the Prognosis and Treatment of Clear Cell Renal Cell Carcinoma

Abstract

Background: Previous studies have shown that RNA N6-methyladenosine (m⁶A) plays an important role in the construction of the tumor microenvironment (TME). However, how m⁶A plays a role in the TME of clear cell renal cell carcinoma remains unclear.

Methods: Based on 23 m⁶A modulators, we applied consensus cluster analysis to explore the different m⁶A modification profiles of ccRCC. The CIBERSORT method was employed to reveal the correlation between TME immune cell infiltration and different m⁶A modification patterns. A m⁶A score was constructed using a principal component analysis algorithm to assess and quantify the m⁶A modification patterns of individual tumors.

Results: Three distinct m⁶A modification patterns of ccRCC were identified. The characteristics of TME cell infiltration in these three patterns were consistent with immune rejection phenotype, immune inflammation phenotype, and immune desert phenotype. In particular, when m⁶A scores were high, TME was characterized by immune cell infiltration and patient survival was higher (p < 0.05). When m⁶A scores were low, TME was characterized by immunosuppression and patient survival was lower (p < 0.05). The immunotherapy cohort confirmed that patients with higher m⁶A scores had significant therapeutic advantages and clinical benefits.

Conclusions: The m⁶A modification plays an important role in the formation of TME. The m⁶A scoring system allows the identification of m⁶A modification patterns in individual tumors, discriminates the immune infiltrative features of TME, and provides more effective prognostic indicators and treatment strategies for immunotherapy.

Figure 1

The flow chart of our study of identifying hallmark genes and candidate agents.

Figure 2

(a) The mutation frequency of 23 m⁶A regulatory factors in 336 patients with clear cell renal cell carcinoma from the TCGA-STAD cohort. Each column represents an individual patient. The bar graph above shows TMB, and the numbers on the right indicate the mutation frequency of each regulator. The bar graph on the right shows the proportions of each variant type. The stacked bar chart below shows the conversion rate in each sample. (b) Frequency of CNV changes in m⁶A modulators in the GSE29609 cohort. The height of the column represents the frequency of change. Delete frequency, blue dot; zoom in frequency, red dot. (c) The expression of 23 m⁶A expression factors between normal tissues and clear cell renal cell carcinoma tissues. Tumor, red; normal, blue. The upper and lower ends of the box represent the interquartile range of values. The line in the box represents the median value, and the red or blue dots represent outliers. The asterisk represents the statistical p value (^∗p < 0.05; ^∗∗p < 0.01; and ^∗∗∗p < 0.001) (d) The position of the CNV change in the m⁶A regulatory factor on the 23 chromosomes of the GSE29609 cohort.

Figure 3

The patterns of m⁶A methylation modification and the biological characteristics of each pattern. (a–d) Using unsupervised cluster analysis to show that ccRCC can be divided into three different genotypes (k = 3). (e) Unsupervised clustering of 23 m⁶A regulatory factors in the clear cell renal cell carcinoma cohort. Survival status, clinical stage, age, project, and m⁶A cluster are used as patient annotations. Red represents the high expression of regulatory factors, and blue represents low expression. (f) Survival analysis based on three m⁶A modification patterns of 555 gastric cancer patients from 1 GEO cohort, including 251 cases in m⁶A cluster-A, 205 cases in m⁶A cluster-B, and 99 cases in m⁶A cluster-C. The Kaplan-Meier curve with a log-rank p-value of 0.011 shows a significant difference in survival between the three m⁶A modification modes. The overall survival rate of m⁶A cluster-A is significantly better than the other two m⁶A clusters. (g) Interaction between m⁶A regulatory factors in clear cell renal cell carcinoma. The size of the circle represents the influence of each adjusting factor on the prognosis, and the numerical ranges calculated by the log-rank test are p < 0.001, p < 0.01, p < 0.05, and p < 0.1, respectively. The purple dots in the circle are prognostic risk factors; the green dots in the circle are prognostic protective factors. The lines connecting the regulators show their interaction, and the thickness shows the relative strength between the regulators. Negative correlations are marked in blue, and positive correlations are marked in red. Writers, erasers, and readers are marked in gray, red, and orange. (h–j) Correlation between the regulators.

Figure 4

(a-b) GSVA enrichment analysis shows the activation status of biological pathways in different m⁶A modification modes. Heatmaps are used to visualize these biological processes. Red represents activated pathways, and blue represents inhibited pathways. The clear cell renal cell carcinoma cohort was used as sample annotation. A m⁶A cluster A and m⁶A cluster B; B m⁶A cluster B and m⁶A cluster C. (c) The abundance of each TME-infiltrating cell in the three m⁶A modification modes. The upper and lower ends of the box represent the interquartile range of values. The lines in the boxes represent the median value, and the black dots represent the outliers. The asterisk represents the statistical p value (^∗p < 0.05; ^∗∗p < 0.01; and ^∗∗∗p < 0.001). (d) The principal component analysis of the transcriptome profile of the three m⁶A modification patterns shows significant differences in the transcriptome among different modification patterns.

Figure 5

(a-b) Function annotation, m⁶A-related genes using GO enrichment analysis, and KEGG enrichment analysis. The color depth of the bubble chart represents the number of enriched genes, and the size of the bubble chart represents the proportion of gene expression. Unsupervised clustering of overlapping m⁶A phenotype-related genes in the C clear cell renal cell carcinoma cohort to classify patients into different genomic subtypes, respectively, is called m⁶A gene clusters A-C. Survival status, clinical stage, age, m⁶A clusters, and gene clusters were used as patient annotations. (c–f) Using unsupervised cluster analysis to show that ccRCC can be divided into three different genotypes (k = 3). (g) A total of 1,152 DEGs were obtained from the three types. (h) Three different m⁶A modified genome phenotypes were discovered through an unsupervised clustering algorithm, and these three clusters were named m⁶A gene clusters A, B, and C. (i) Kaplan-Meier curve showed that the m⁶A modified genome phenotype was significantly correlated with the overall survival of 555 patients in the clear cell renal cell carcinoma cohort, including 296 gene cluster A 120 gene cluster B and 139 gene cluster C (p < 0.0001, log-rank test). (j) The expression of 23 m⁶A regulatory factors in 3 gene clusters. The upper and lower ends of the box represent the interquartile range of values. The line in the box represents the median value, and the red or yellow dots represent the outliers. The asterisk represents the statistical p value (^∗p < 0.05; ^∗∗p < 0.01; and ^∗∗∗p < 0.001). One-way analysis of variance is used to test the statistical differences between the three gene clusters.

Figure 6

(a-b) Waterfall plots of tumor somatic mutations established by those with high m⁶A score (a) and low m⁶A score (b). Each column represents a single patient. The bar graph above shows TMB, and the numbers on the right indicate the mutation frequency of each gene. The bar graph on the right shows the proportion of each variant type. (c) The Sankey diagram shows the changes in m⁶A clusters, gene clusters, m⁶A score, and survival status. (d) Survival analysis of low (62 cases) and high (493 cases) m⁶A score patient groups in the TCGA-KIRC cohort, using Kaplan-Meier curve (HR, 1.81 (1.26–2.62); p < 0.001, log-rank test). (e) Spearman's correlation analysis was used to analyze the correlation between m⁶A score and immune cells in the clear cell renal cell carcinoma cohort. Negative correlations are marked in blue, and positive correlations are marked in red. (f) Using the Kaplan-Meier curve to analyze the survival rate of patients with low (254 cases) and high (74 cases) tumor mutations in the clear cell renal cell carcinoma cohort (HR, 1.81 (1.26–2.62); p < 0.001, log-level test). (g) Kaplan-Meier curve is used to analyze survival by m⁶A score and TMB score. H high; L low. (h-i) Differences in m⁶A score between the three gene clusters and the three types in the clear cell renal cell carcinoma cohort. The Kruskal-Wallis test was used to compare the statistical differences between the three gene clusters (p < 0.001).

Figure 7

(a–f) Kaplan-Meier curve is used to show the difference in survival rate between the high group and low group in age and clinical stage (p < 0.05 is meaningful). (g–i) The difference in the expression levels of the three immune markers in the high- and low-score groups (p < 0.1 is meaningful). (j) The distribution of scoring scores between the high and low groups in ccRCC, and the difference between the high and low groups, and the black bars represent the median. (k–n) No anti-PD-1/L1 immunotherapy and anti-CTLa4 immunotherapy, anti-PD-1/L1 immunotherapy alone, anti-CTLa4 immunotherapy alone, anti-PD-1/L1 immunotherapy and anti-PD-1/L1 immunotherapy, and anti-PD-1/L1 immunotherapy alone were separately received. The different efficacy of CTLa4 immunotherapy in the high- and low-score groups when the two combined immunotherapy is combined, and the black bars represent the median.