The mRNA transcriptomic profile and corresponding clinical data of ccRCC samples were downloaded from The Cancer Genome Atlas (TCGA), including 72 normal controls and 539 cancer samples. The detailed information of patients’ information has been showed in the Table 1. Duplicate samples and those with an overall survival time of fewer than 30 days were excluded. Later, when building the LASSO regression model, we retained TCGA samples containing both RNA-seq profiles and survival data for analysis. We installed a 3:2 ratio to randomly group these samples as the training cohort (n = 315) and internal validation cohort (n = 208). Furthermore, expression profiling and clinical information on the validation cohorts were downloaded from the ICGC and E-MTAB databases. All data were processed through R (version 4.0.3) and converted to FPKM value format. In addition, we obtained eight ccRCC datasets from the GEO database, transformed the data to the Z score, and validated the expression levels of KNCJ16. We built an outline diagram for the whole process of the study (Figure 1).

We first screened and documented ion channel-related genes (ICRGs) through an extensive literature review. The detailed list of the ICRGs is given in Supplementary Table S1. The differentially expressed genes (DEGs) between tumor and normal samples were screened using R studio’s “Limma” package. |Log2FC| > 1, p < 0.05, and false discovery rate (FDR) < 0.05 were set as the cutoffs for the DEGs in both TCGA and ICGC cohorts. An interaction network between the DEGs was mapped by GeneMANIA website (Warde-Farley et al., 2010). The DEGs that were significant in both cohorts were used for subsequent analysis. Univariate Cox regression was used to identify the prognostic ICRGs. Genes with a p-value less than .05 in Cox regression were identified as prognostic genes for further LASSO regression analysis.

Based on the previously obtained DEGs, we used both the least absolute shrinkage and selection operator (LASSO) and support vector machine recursive feature elimination (SVM-REF) to identify diagnostic models in ccRCC and normal tissues (Xia et al., 2020). The genes obtained by the two methods were taken as intersections and validated with ROC curves. The expression of intersecting differential genes in ccRCC and normal tissues was also explored in the Human Protein Atlas (HPA) database.

To subtype patients and thus proceed with individualized treatment, we performed molecular subtyping of patients in TCGA based on prognostic ICRGs using the ConsensusClusterPlus package (Wilkerson and Hayes 2010). The RNA matrix was calculated using the k-means algorithm on a fraction of probes (0.8) with 1,000 repeats. Meanwhile, unsupervised clustering and corresponding representative figures were generated with the ggplot2 package (Ito and Murphy 2013). Then, a single-sample Gene Set Enrichment Analysis (ssGSEA) was performed to identify immune cell infiltration between the two clusters. The potential molecular enrichment of the two clusters was annotated with the ClueGO plugin (Bindea et al., 2009). The pRRophetic algorithm was performed to calculate the half-maximal inhibitory concentration (IC50) of chemotherapeutic drugs in clusters (Geeleher et al., 2014). At the same time, we applied The Cancer Immunome Atlas (TCIA) database to observe the response to immunotherapy between the two clusters.

The pathway enrichment analyses were conducted using the ClusterProfiler packages to explore the potential functions, and the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were chosen as the background databases (Yu et al., 2012). The Gene Ontology (GO) analysis identified biological processes (BPs), molecular functions (MFs), and cell components (CCs) enriched for each set of genes.

To assess the value of ICRGs as a guide to the prognosis of ccRCC, we first divided the samples in TCGA into a training and validation set on a 3:2 basis while employing the E-MTAB-1980 database as an external validation set. Furthermore, we trained the prognostic ICRGs on LASSO with maxit = 1,000 and constructed a 10-gene prognostic model. ROC curves and K–M survival plots were condemned to assess the predictive value of prognostic models using survival and time ROC packages.

To better assess clinical prognosis and provide timely intervention, we combined risk scores and clinical parameters to construct the nomogram. The “rms” package was used to establish the nomogram of multivariable models for 1-, 3-, and 5-year OS prediction. Furthermore, the calibration curve was carried out for the constructed nomogram, and the discriminative performance is quantified using the area under the curve (AUC).

In the results of this study, continuous variables are described by mean ± S.D, while categorical variables are characterized by frequency (n) and proportion (%). Univariate and multivariate Cox regression analyses were used to detect independent factors. The survival analysis evaluated the correlation between the characteristics and overall survival by LASSO regression models. Kaplan–Meier survival curves and time-dependent ROC curves were drawn and compared among subgroups. In addition, statistical analyses were performed using R version 4.0.3 and GraphPad 9.0. A p-value of less than .05 was considered statistically significant.

The expression levels of ICRGs were analyzed between ccRCC and normal samples in both TCGA and ICGC databases with a p-value < 0.05 and absolute log₂FoldChange >1. The volcano plot exhibited the DEGs in TCGA and identified the five most significantly up- and downregulated genes (Figure 2A). Afterward, we considered 113 genes altered considerably in both datasets as DEGs in ccRCC (Figure 2B). Subsequently, to understand the correlation between these differential genes, we constructed a gene interaction network by GENEMANIA (Figure 2C). The interaction network demonstrated a variety of potential connections between DEGs, including potential interaction links such as co-expression and co-localization. Then, the GO enrichment analysis showed an enrichment of GO terms indicative of functional maturation, such as ion transport or ion channel composition, which is consistent with the function of these ICGRs (Figure 2D). In addition, the analysis based on KEGG pathways identified several significant enrichments, such as nicotine addiction and the calcium signaling pathway (Figure 2E).

The DEGs were obtained and their potential associations and enriched pathways were analyzed. Next, we explored whether ICRGs could be used as a diagnostic model to determine ccRCC. We screened 10 and nine ICRGs by LASSO regression and SVM-REF analyses, respectively (Figures 3A, B). We then selected the intersecting genes of the two methods as diagnostically relevant genes, which were KCNJ1, KCNJ10, KCNJ16, TRPV6, and GABRA2 (Figure 3C). Subsequently, we evaluated the predictive diagnostic value of these five ICRGs by AUC values in TCGA and ICGC databases, and the results showed a good predictive value in ccRCC or RCC (Figures 3D, E). Subsequently, we explored the protein expression levels of these genes in the HPA database. It could be seen that KCNJ1, KCNJ10, KCNJ16, and TPRV6 expression levels were all significantly upregulated in normal tissues compared with tumor tissues and distributed in different normal tissue structures (Figure 3F). Due to the lack of relevant IHC staining for GABRA2 in the database, the IHC image of GABRA2 has not been exhibited. Corresponding informations has been changed in the manuscript too.

Molecular subtyping could provide insights into tumor biology for potential therapeutic targets. So, in this section, we explored whether ICRGs could be used in molecular subtyping for ccRCC patients. First, by performing the Univariate Cox analysis, we screened out 33 ICRGs significantly associated with OS and defined them as prognostic ICGRs. The mRNA expression level of 33 prognostic ICRGs between ccRCC and normal samples was exhibited (Figure 4A). We then conducted the ConsensusClusterPlus package and identified two clusters (C1 and C2), according to the consensus distribution function plot (Figures 4B, C). Principal component analysis (PCA) showed good separation between the two classifications and the survival showed a worse prognosis for patients in C2 than in C1. We then drew a heatmap of the composite clinical information based on the clusters. The gene expression differences between the two clusters and the corresponding clinical grading and staging distribution can be apparently shown (Figure 4E).

To explore potential features between the two clusters, we identified gene expression differences between the clusters by using the limma package and performed functional enrichment using the ClueGO plugin. We found that patients in cluster one were significantly associated with connexin complex, nicotine addiction, and positive regulation of potassium ion transmembrane transport, while patients in cluster two were significantly correlated with TRP channels (Figure 5A). Subsequently, we explored the sensitivity of the two clusters to drugs commonly used in clinically advanced ccRCC by pRRophetic. The results showed that the patients in cluster1 were more sensitive to Axitinib and Cytarabine, while the patients in cluster 2 were more sensitive to the remaining four drugs (Figure 5B). Subsequently, differences in immune cell infiltration were explored between our two clusters. Immune cell infiltration was more abundant in cluster 2 (Figure 5C). Also, we explored the responsiveness to immunotherapy in the TCIA database. More patients in cluster two were PD-1 positive or PD-1/CTLA-4 double positive and may have benefited more from immunotherapy (Figure 5D).

Based on the aforementioned prognostic ICRGs, we further evaluated the prognostic value of ICRGs in ccRCC patients. The model coefficients were solved by LASSO regression (Figure 6A, B). The detailed information on the ICRGs from the signature is shown in Table 2. The risk score is constructed according to the following formula mentioned in the Methods section. Based on the calculated median risk score cut-off, patients were divided into high- and low-risk groups. Furthermore, the Kaplan–Meier log-rank test and the time-dependent ROC curve were used to evaluate the predictive ability and accuracy of the prognostic signature. The outcome of the Kaplan–Meier log-rank test showed that the high-risk group had a significantly worse OS compared with the low-risk group in TCGA training set (Figure 6C), TCGA validation set (Supplementary Figure S1A), and E-MTAB validation set (Supplementary Figure S1C). Meanwhile, the time-dependent ROC curve proved the 1-, 3-, and 5-year predictive accuracy of the signature for OS (Figure 6D; Supplementary Figures S1B, D). In addition, the risk score distribution, survival status, and expression of ICRGs from the signatures are exhibited in TCGA training set, TCGA validation set, and E-MTAB validation set (Figure 6E; Supplementary Figures S1E, F).

To better predict the prognosis of patients, we combined clinical factors and risk signatures to construct the nomogram. Univariate Cox regression analysis was performed to screen out the significant factors with OS (Figure 7A). Afterward, multivariate Cox regression analysis showed that the AJCC stage, age, and risk score could be the independent factors for OS of ccRCC patients (Table 3). Combining the risk score and clinical parameters, we established the nomogram with a C-index of 0.774 (Figure 7B). Then, we performed the calibration curves to verify the predictive efficacy of the nomogram for 1-, 3-, and 5-year OS (Figure 7C). We also performed ROC curves, and the AUC values indicated that the predictive efficacy of the nomogram was significantly better than that of the present AJCC stage or ISUP grade and the prognostic signature alone (Figure 7D).

We then analyzed the KNCJ16 gene, an ICRG that we had previously derived to be significant in both diagnostic and predictive models. First, we found that the KCNJ16 expression level was significantly altered in a variety of cancers by pan-cancer analysis, using the SangerBox and was more apparent in ccRCC (Figure 8A). We then compared the expression levels of KCNJ16 in cancer and normal tissues in eight GEO databases and found that the expression levels of KNCJ16 were significantly upregulated in normal tissues (Figure 8B). Based on the aforementioned results, we wanted to clarify exactly in which structure of the kidney that the KCNJ16 expression levels are upregulated, and we explored this with single-cell data from the HPA database. We found that the KCNJ16 gene was mainly upregulated in the proximal tubular cells and collecting duct cells (Figure 8C). Meanwhile, we screened out the relationship between KNCJ16 and immune cell infiltration using the TIMER database. Consistent with the HPA database, the correlation between KCNJ16, T cells, and B cells was not significant, while mainly affecting macrophage components (Figure 8D). Furthermore, we explored the potential biological functions of KCNJ16. GSEA enrichment analyses showed that high KCNJ16 expression was significantly associated with protein secretion, p53 pathway, and so on (Figure 8E). Meanwhile, GO and KEGG enrichment analyses indicated that KCNJ16 was closely related to the transmembrane transport function and cell structure maintenance (Figure 8F).