A total of 28 RNA-seq transcriptional profiling of normal samples were downloaded from the GTEx (Genotype-Tissue Expression) dataset (https://gtexportal.org/). A total of 1021 RNA-seq transcriptome profiling (128 normal samples and 893 RCC samples), 872 DNA methylation data (205 normal samples and 667 tumor samples), and corresponding clinical information of RCC were downloaded from TCGA (The Cancer Genome Atlas) dataset (https://gdc.cancer.gov/).

We standardized RNA-seq transcriptional profiling by using the “limma” R package, and the Wilcoxon rank-sum test was utilized to identify DEGs (Ritchie et al., 2015; Zhang et al., 2020). The false discovery rate (FDR) < 0.05 and |log₂ fold change (FC)| > 2 were taken advantage of as cutoff criteria. The “pheatmap” R package was used to plot the heatmaps (Li et al., 2020).

Gene expression data and DNA methylation data were standardized by using the “limma” R package (Ritchie et al., 2015). DNA methylation-driven genes (MDGs) were identified using the “methylMix” R package (Gevaert, 2015; Zhang et al., 2020). DNA methylation data and paired gene expression data were integrated and analyzed jointly to identify the DNA methylation status negatively correlated with the gene expression of a particular gene, indicating that the gene is a DNA methylation-driven gene (Cedoz et al., 2018). Inclusion criteria were set to the correlation of methylation data and corresponding gene expression data of DEGs less than −0.3, |log₂FC| > 0 and adjust p < 0.05. The differentially expressed DNA methylation-driven genes (DEMDGs) were obtained by intersecting MDGs and DEGs for further analysis (Zhang et al., 2020).

The PPI network was established by the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/) and visualized by Cytoscape software (v3.8.2) (Shannon et al., 2003). To generate an interaction network, 146 DEMDGs were uploaded to the STRING database. The obtained networks were downloaded in a tabular format and uploaded to Cytoscape for network visualization.

The immune cell infiltration levels in the samples were evaluated by the “CIBERSORT” R package, and the inclusion criteria were p < 0.05 (Newman et al., 2015). The stromal score and immune score were also estimated by the “estimate” R package (Jia et al., 2018). Boxplots were plotted using the “ggpubr” R package (Whitehead et al., 2019).

We capitalized on the k-means clustering algorithm in the “ConsensusClusterPlus” R packet to perform clustering (Wilkerson and Hayes, 2010). Here, we performed the partition around medoids clustering algorithm and Euclidean distance. The cluster number was tested from 2 to 9, and the optimal one was selected to produce the most stable consensus matrix and the least overlapping cluster assignments across permuted clustering runs. We implemented it with the R package ConsensusClusterPlus.

Utilizing the “glmnet” R package, “survival” R package, and “survminer” R package, the prognostic model was constructed (Armbruster et al., 2019; Wang et al., 2021; Xi et al., 2021). The risk scores were calculated according to a linear combination of the gene expression levels weighted by the regression coefficients from the multivariate Cox regression analysis (Wei et al., 2020). The “survivalROC” R package was utilized to validate the stability of the prognostic model (Huang et al., 2017). The Kaplan–Meier survival curves were carried out to assess the survival time between high- and low-risk score RCC patients (Lawrie et al., 1982). We validated the accuracy of the optimal cut-off value by the principal component analysis (Kim et al., 2021). The Human Protein Atlas database (https://www.proteinatlas.org/) was explored to investigate the expression level of the prognostic genes.

To evaluate the expression level of mRNA, total RNA was extracted from the cell using TRIzol (Invitrogen). One microgram of total RNA was used as a template for cDNA synthesis using a PrimeScript RT Reagent Kit (Takara). Real-time quantitative PCR (qRT-PCR) was performed in a reaction mixture containing SYBR Green (Takara) with the CFX96 Touch real-time PCR detection system (Bio-Rad). The expression level of related mRNAs was calculated using 2^−ΔΔCT, and the related GAPDH mRNA expression was used as an endogenous control. The primer sequences involved in this study are shown in Supplementary Table S1. Each PCR reaction was performed in triplicates.

We used the “survival” R package to evaluate the independence of the prognostic model from clinical features via univariate and multivariate Cox regression analyses (Qi et al., 2021). The significant levels were set to p < 0.05, and hazard ratios (HRs) with 95% CIs were also calculated.

The nomogram was constructed utilizing the “rms” R package (Pond et al., 2014). The web-based survival rate calculator was established using the “shiny” and “DynNom” R packages to predict cancer-specific survival rates dynamically (Sun et al., 2020; Bakin et al., 2021). Calibration curves, which plot the average Kaplan–Meier evaluation according to the corresponding nomogram for 1-, 3-, or 5-year predicted overall survival, are provided to estimate the accuracy of the nomogram.

Gene set enrichment analysis (GSEA) was performed using GSEA4.0 (https://www.gsea-msigdb.org/gsea/index.jsp/). The annotated gene set files (c2.cp.kegg.v7.4. symbols.gmt gene set and c5.go.bp.v7.4. symbols.gmt) were considered as the reference gene set. The inclusion criteria were p < 0.05 and FDR < 0.25. The column diagrams were plotted by GraphPad Prism 7 (Elzein et al., 2021).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed using the “enrichplot” R package, “org.Hs.eg.db” R package, and “clusterProfiler” R package (Yu et al., 2012; Zhang et al., 2019). The inclusion criteria were set to p < 0.05 and q < 1. The results were visualized by the “ggplot2” R package (Sun et al., 2011).

All statistical tests were performed by R statistical software (version 4.0.3) (http://www.r-project.org/) using Mann-Whitney testing for continuous data and Fisher’s exact testing for categorical data. The correlation between two continuous variables was measured by Pearson’s correlation coefficient. The hazard ratio (HR) and 95% confidence intervals (CI) were estimated by a Cox regression model using the survival package. Survival analysis was carried out using Kaplan–Meier methods. Differential methylation was calculated from mean (β-value-cancer)—mean (β-value- normal). The differences in variables among different groups were compared by means of the Student’s t-test. p < 0.05 was considered to be statistically significant.

The research process of the study was shown in Figure 1. Based on the Wilcoxon rank test, a total of 5,198 DEGs were selected from 28 RNA-seq transcriptome profiling of normal samples in the GTEx dataset and 1021 RNA-seq transcriptome profiling (893 RCC samples and 128 normal samples) in TCGA dataset (FDR <0.05 and |log₂FC| > 2). The heatmap shows the expression of DEGs between RCC samples and normal samples (Figure 2A). We screened the 270 methylation-driven genes (MDGs), whose methylation status negatively correlated with expression levels (Cor < −0.3, |log₂FC| > 0 and adjust p < 0.05) (Supplementary Table S2). The heatmap shows the expression of MDGs between RCC samples and normal samples (Figure 2B). Then, 270 MDGs and 5,198 DEGs were intersected to obtain 146 DEMDGs for further analysis (Figure 2C). We further visualized the methylation levels (Figure 2D) and gene expression levels (Figure 2E) of 146 DEMDGs in RCC samples and normal samples. The comprehensive landscape of DEMDG interactions and DEMDG connection for RCC patients was depicted with the DEMDG network (Figure 2F). The aforementioned results found 146 DEMDGs in RCC and normal samples. These abnormal DEMDGs were interconnected and may be involved in the occurrence and development of RCC.

To select the optimized cluster number, we calculated the k-means clustering algorithm with the ConsensusClusterPlus R packet. K = 2 was identified with the optimal clustering stability (Figures 3A–D). Then, we analyzed the methylation levels and gene expression levels of 146 DEMDGs, as well as the clinical features of paired patients. There were significant differences in the methylation levels and gene expression levels between cluster 1 and cluster 2, and the clinical features were evenly distributed in two clusters (Figures 3E,F). The RCC patients in cluster 2 (n = 419) had better overall survival (OS) than the patients in cluster 1 (n = 435, p < 0.001) (Figure 3G).

Immune checkpoint inhibitors (ICIs) are administered for the treatment of RCC. We investigated whether the two clusters were related to ICI-related biomarkers. The results showed that cluster 1 was positively correlated with the high expression of LAG3 (p < 0.001), CD160 (p < 0.001), HAVCR2 (p < 0.001), CTLA4 (p < 0.001), and TIGIT (p < 0.001), and the stromal score and immune score were significantly higher in cluster 1 than in cluster 2 (p < 0.001) (Figure 4A). The abundance of naive B cells (p < 0.001), CD8 T cells (p < 0.001), CD4 memory-activated T cells (p < 0.001), follicular helper T cells (p < 0.001), gamma delta T cells (p < 0.001), and macrophages M1 (p < 0.001) was significantly higher in cluster 1 than in cluster 2 (Figure 4B). The higher immune infiltration level corresponded to cluster 1, and the lower immune infiltration level corresponded to cluster 2 (Figure 4C). The RCC patients in high-immune score groups had a worse OS than the patients in low-immune score groups (p < 0.001) (Figure 4D). Furthermore, we accessed the correlation of immune cells in cluster 1 and cluster 2. In cluster 1, the positive correlation between CD8 T cells and follicular helper T cells was the strongest, in which the correlation coefficient was 0.55. The correlation coefficient between CD8 T cells and CD4 memory-resting T cells was −0.66, which was the lowest negative correlation (Figure 4E). However, in cluster 2, the cells with the strongest negative correlation were activated CD8 T cells and macrophages M2, in which the correlation coefficient was −0.46 (Figure 4F). These results showed that the two clusters based on 146 DEMDGs were closely associated with prognosis and immune status in RCC patients.

To explore the possible reasons causing worse OS in cluster 1, we selected 106 DEGs from two clusters (Cor < −0.3 and |log₂FC| > 1). We used a heatmap visualizing the gene expression levels of 106 DEGs in RCC and normal samples (Supplementary Figure S1). We performed GO and KEGG analysis to analyze underlying functions and pathways of 106 DEGs (p < 0.05). Results of GO analysis were significantly enriched in regulation of T-cell activation, T-cell proliferation, positive regulation of leukocyte proliferation, positive regulation of T-cell proliferation, positive regulation of cell–cell adhesion, etc. (Figure 4G; Supplementary Figures S2A, S2B). Results of the KEGG pathways were significantly enriched in pathogenic Escherichia coli infection, natural killer cell-mediated cytotoxicity, viral myocarditis, one carbon pool by folate, cytokine–cytokine receptor interaction, etc. (Figure 4H). The aforementioned results indicated that the 106 DEGs were in close contact with the immune microenvironment, which may be the cause for the OS difference between the two clusters.

To determine the prognostic value of 146 DEMDGs, univariate Cox regression analysis, LASSO, and multivariate Cox regression analysis were used to identify them. Subsequently, the prognostic model was constructed based on eight independent and prognostic DEMDGs (including LCP2, PPP1R18, APOL1, FMNL1, CLDN7, NMI, FAXDC2, and SHC1) (Figures 5A–D). We analyzed the association between the gene expression and the survival of the patients. The patients’ OS with high expressions of LCP2, PPP1R18, APOL1, FMNL1, NMI, and SHC1 were worse than the low-expression groups (p < 0.05), while the patients’ OS with high expressions of CLDN7 and FAXDC2 were better than the low-expression groups (p < 0.05) (Supplementary Figure S3). The methylation levels of eight DEMDGs were inversely correlated with their expression levels (p < 0.001) (Supplementary Figure S4). The risk score was calculated as follows = LCP2 * (−3.619) + CLDN7 * (−1.577) + FAXDC2*(-0.977) + APOL1*(1.365) + NMI *(1.686) + PPP1R18 * (1.695) + SHC1 * (2.007) + FMNL1 * (2.313). The coefficients of each gene are shown in Table 1.

RCC patients were split into high- and low-risk groups according to the optimal cut-off value of the risk score (cutoff = 1.78) (Figure 5E). The AUC of the ROC curves was 0.738, 0.673, and 0.703 within 1, 3, and 5 years, which demonstrated that the risk score had a good prognostic value (Figure 5F). The distributions of the risk score in high- and low-risk groups were shown in Figure 5G. Patients’ mortality risk increased with increasing risk scores (Figure 5H). The survival curve was carried out to assess the survival time between high- and low-risk score groups. The survival time of high-risk groups was significantly worse than the low-risk groups (p < 0.001) (Figure 5I). RCC samples were clearly structured in two different groups by the principal component analysis, which suggested our study could significantly reflect the prognosis differences of RCC patients (Figure 5J).

To verify the expression levels of eight DEMDGs in RCC cells, we used RT-qPCR analysis to detect human embryonic kidney-293 (HEK-293) and human RCC cell lines (786-O and OS-RC2 cells) (Figures 6A–H). Among them, six DEMDGs (LCP2, PPP1R18, APOL1, FMNL1, NMI, and SHC1) were upregulated in both RCC cells, combined with Supplementary Figure S3; their high expression was associated with poor survival, considering that they might play a role as proto-oncogenes. FAXDC2 was downregulated in both RCC cells. However, CLDN7 was downregulated in 786-O cells and upregulated in OS-RC2 cells. To determine the clinical relevance of these eight gene expressions, HPA clinical specimens were used to analyze the proteins’ expression encoded by these eight genes (Figures 6I–O). Relative to its expression level in normal kidney tissue, SHC1 were strongly positive, while LCP2, PPP1R18, and NMI were moderately positive in RCC tissues. FMNL1 and CLDN7 were not detected in RCC tissues. APOL1 was expressed more abundantly in the normal tissue than in malignant. FAXDC2 was not found on the website. However, our results on APOL1 are contrary to the database. Our team speculated that the main reason was that the data from TCGA database came from all pathological types of RCC, hence the inconsistent results.

To further annotate functions enriched in the high- and low-risk groups, GSEA was queried to confirm the signaling pathways in which the genes are enriched. The results are represented in Figures 7A–D. The following biological processes were enriched in the high-risk groups: collagen fibril organization, lymphocyte activation, negative regulation of B-cell activation, regulation of the leukocyte apoptotic process, and regulation of leukocyte proliferation. The following signaling pathways were enriched in the high-risk groups: cytokine–cytokine receptor interaction. To clarify the possible molecular mechanism of eight prognosis-related DEMDGs, we also performed GO and KEGG pathway analysis. Results of the GO analysis were significantly enriched in the sterol metabolic process, Fc-epsilon receptor signaling pathway, Fc receptor signaling pathway, positive regulation of megakaryocyte differentiation, interleukin-2 mediated signaling pathway, etc. (Figure 7E). Results of KEGG analysis were significantly enriched in the natural killer cell-mediated cytotoxicity, African trypanosomiasis, Fc epsilon RI signaling pathway, prolactin signaling pathway, and chronic myeloid leukemia (Figure 7F). These results suggested that RCC patients’ prognosis might be impacted by the above biological functions and signaling pathways.

To determine whether the risk score could be presented as an independent prognostic factor for RCC patients, we employed univariate and multivariate Cox proportional hazard regression analyses. In the univariate analysis and multivariate analysis, the risk score, age, and M stage showed pronounced effects on the RCC prognosis (p < 0.05) (Figures 8A, B). Furthermore, we constructed a nomogram with the significant variables in the multivariate analysis (Figure 8C). Results suggested that the risk score had a significant influence on survival prediction. The 1-, 3-, and 5-year predicted calibration curves also respectively suggested that the model had a good prediction accuracy (Figures 8D–F). We also established a dynamic web-based survival rate calculator (http://127.0.0.1:3496), which could individually predict the survival of patients according to their clinical features and risk score. For example, the 3-year cancer-specific survival rate was approximately 76% (95% CI 42–72%) for patients with low risk, M0 stage, and aged <65 years (Figures 8G, H).

To further analyze the relationship between the prognostic model and tumor-infiltrating immune cells, we performed a detailed Spearman correlation analysis, and the result was presented with the lollipop shape (Figure 9A). High-risk groups were more positively correlated with tumor-infiltrating immune cells, including CD8⁺ T cells, macrophage M1, B cells, monocytes, myeloid dendritic cells, regulatory T cells, and myeloid dendritic cells. The detailed results are shown in Supplementary Table S3. We also attempted to identify associations between the prognostic model and the efficacy of six common chemotherapeutic drugs for the treatment of RCC. The high-risk score was associated with the lower half-maximal inhibitory concentration (IC50) of chemotherapeutics such as temsirolimus (p = 5.4e−15), sunitinib (p < 2.22e-16), pazopanib (p = 0.011), sorafenib (p = 0.015), bleomycin (p = 0.00017), and axitinib (p = 2.2e−13) (Figures 9B–G). The aforementioned results showed that this model closely correlated with tumor-infiltrating immune cells and drug susceptibility.