The ccRCC SCQ dataset GSE152938 was downloaded from the GEO database, including 1 normal kidney sample, 2 ccRCC samples, 1 chromophobe renal cell carcinoma sample and 1 papillary renal cell carcinoma sample. Because this article was designed to study the prognosis of patients with ccRCC, we removed other types of samples. Next, we perform quality control on the SCQ data, we selected genes expressed in at least three cells and cells with total gene expression between 300 and 3000 for the next analysis. And, cells with mitochondrial gene expression greater than 5% of total gene expression were also excluded. The transcriptome RNA-seq data and its corresponding clinical information were acquired from the TCGA database, comprising 539 ccRCC data and 72 normal data. And, to ensure the accuracy of the study, 530 ccRCC samples that included complete clinical information were selected for further analysis.

First, we normalized the ccRCC SCQ data filtered in the previous step by the method of “LogNormalize”. Then, due to the sheer volume of cells, we classified them by the marker genes expressed by each cell, and merged the similar categories, through the method of principal component analysis (PCA) dimension reduction. Finally, with the help of the function of “SingleR”, we annotated cell types according to their marker genes, and we used the “FindAllMarkers” function to obtain marker genes for different cell types.

WGCNA is a systematic statistical approach that can group genes with analogous expression patterns and illustrate the relationship between genes of a particular group and specific traits (19). In our study, we used this method to obtain the set of genes associated with clinical traits. First, we performed an initial screening of the samples, excluding non-renal cancer patients and genes with small fluctuations. Then, to improve the accuracy of the screening, we transformed the adjacency matrix into a topological overlap matrix (TOM) and set the minimum group size to 30. Finally, we merge similar groups and output the resulting graph and data.

Firstly, we combined the above-obtained genes with TCGA transcriptome data to obtain the immune subtype-related gene expression data. Subsequently, we merged the expression data with their survival status and performed a univariate Cox analysis to identify genes associated with prognosis. Then, the genes were further selected by the way of Lasso regression analysis, and through this we can get the model genes. Finally, we calculated the risk score for each ccRCC patient based on the formula and with the help of the median risk score, we were able to classify the patients into two risk groups.

We analyzed whether the risk score was an independent prognostic factor by the measure of Cox analysis. Then, we assessed the predictive effect of the model by plotting the survival curves of the ccRCC patients, and we make the most of the ROC curve to evaluate the accuracy and sensitivity of this model. In the last, we plot the patient’s survival status on an axis with the risk score as the horizontal coordinate to give a better visualisation of each patient’s survival status.

With the help of the results of 7 kinds of immune infiltration in ccRCC downloaded from the TIMER database, we showed the difference in the level of immune infiltration between ccRCC patients in two groups. Meanwhile, we also investigated the difference in the expression level of immune checkpoint-related genes between ccRCC patients in two groups.

In our study, by combining each patient’s risk score of our model with clinical information, we successfully constructed a nomogram that can predict the patient’s risk of death. Then we used ROC to assess the accuracy of nomogram in predicting patient outcomes

CcRCC cell-lines 786-O was purchased from the Chinese Academy of Sciences Committee on Type Culture Collection Cell Bank (Shanghai, China) and were cultured in RPMI-1640 medium (Gibco, USA) supplemented with 10% FBS (Gibco, USA). siRNA VMP1 and siRNA negative control were purchased from RiboBio (Guangzhou, China), and transfected with Lipofectamine 2000 reagent (Invitrogen, CA, USA).

Total RNAs of 786-O were extracted using the TRIzol reagent. Then, we made use of a reverse transcription kit from (vazyme, China) to obtain cDNA. We detected the relative expression level of the target gene by the measure of qRT-PCR based on the 2-ΔΔCt method.

5-ethynyl-29-deoxyuridine (EdU) assay was performed based on the manufacturer’s instructions (RiboBio, Guangzhou, China). The 786-O cells were first inoculated in 24-well plates, followed by incubation with EdU reagent for 2h. Finally, after labelling the DNA with 2-(4-Amidinophenyl)-6-indolecarbamidine dihydrochloride (DAPI), the cell images were inspected under a fluorescent microscope.

To assess the migratory capacity of the cells, we adjusted the 24-well plates (Nset, China) by transwell culture chambers (Corning, USA). Cells were inoculated into 200μL of the medium in the upper chambers without serum. The lower layer of the chamber is 700μL of medium containing 10%FBS. For cell invasion ability, pre-lay a layer of Matrigel over the chambers, the rest of the steps are the same as above. After 24h incubation in the cell incubator, the medium was discarded and the cells were wiped from the inside of the bottom of the chambers using a cotton swab. Finally, after fixing them with methanol and staining the cells at the bottom of the chamber with crystal violet, images of the cells were taken using a microscope.

Cells in logarithmic growth phase were inoculated in 6-well plates. When the cells reached about 90%, 200μL tips were used to draw 2 vertical lines along the vertical direction, and after washing out the cell debris, the complete medium was replaced with a medium containing 1% FBS. After 24 hours, the distance the cells migrated to the scratched area was carefully observed under the microscope and this was used to test the migration ability of the cells.

Cells were first digested with trypsin and washed twice with PBS. Then, according to the Annexin V-FITC/PI Apoptosis Detection Kit (vazyme, China) guidelines, cells were incubated with Annexin V-FITC and PI in a dark environment for 10 minutes. Finally, the rate of apoptosis was measured using flow cytometry.

Bioinformatics analysis was conducted using the R software (V. 4.1.2). The quantification and graphing of the experiment data was conducted using Image J software (V.1.8.0) and GraphPad Prism (V.9.0). All measurement data are shown as the mean ± SD. The data differences between the two groups were analyzed by Student’s t-test and P-values less than 0.05 were considered significant in all tests.

As illustrated in Figure 2A , we found a relatively even distribution of cells with gene expression levels between 300 and 3000, and the mitochondrial genes of the majority of cells were <5%. Using the above criteria to screen the cells, we successfully obtained 1759 cells. Figure 2B showed that these cells were uniformly distributed in the ccRCC samples and that the gene expression levels were positively associated with the amount of gene expression (0.88). This indicated that the screened cells are suitable for further analysis. Figure 2C showed the ten most variable genes in selected cells, including JCHAIN, RGS5, ENPP2 and MZB1.

After PCA descending treatment, these cells were divided into 11 clusters. In Figure 3A , we could find 10 highest expressed genes in each cluster. In Figure 3B , we could find the distribution of these 11 clusters. With the help of the function of “SingleR”, we annotated cell types according to their marker genes, and the clusters associated with immune cells are 0, 1, 4, 6, 7, 9 ( Figure 3C ). We then used the “FindAllMarkers” function to acquire 858 immune subtype-related genes.

In TCGA cohort, with the help of WGCNA, we got the gene modules related to the patient’s survival status. By using a soft threshold of 4 and a minimum module gene count of 30, we succeeded in obtaining 3 modules related to clinical traits ( Figures 3D, E ). Because we wanted to analyze patients’ prognosis, we selected genes associated with the patients’ survival status for further analysis.

First, a differential analysis was performed based on genes obtained in the previous step in the TCGA cohort to obtain the differential genes in the tumor and normal groups. Then, as shown in Figure 4A , we succeeded in obtaining 66 genes related to the prognosis of ccRCC patients through univariate Cox analysis, 63 of which had a hazard ratio (HR) > 1. In the last, after randomizing patients into the training and validation set, we carried out Lasso regression analysis on these 64 genes, and the result showed when the number of genes included is 5, the gene contraction tended to be stabilized and the partial likelihood deviation was minimized ( Figures 4B, C ). We finally obtained 5 model genes, including IFI30, CEBPB, VMP1, ATP1B1, and FKBP11, and we found that gene ATP1B1 was highly expressed in normal patients, while IFI30, CEBPB, VMP1, and FKBP11 were highly expressed in ccRCC patients ( Figure 4D ). The names and the coefficients of the prognostic genes were listed in Table 1 . Risk score = IFI30*0.252 + CEBPB *0.050 + VMP1*0.041 + ATP1B1*(-0.007) + FKBP11*0.302. We then used median patient risk values to classify patients into two groups for further analysis.

First, to explore whether risk scores were an independent factor of influence for ccRCC patients, we performed univariate and multivariate Cox regression on age, gender, grade, stage and risk score in ccRCC patients. The presentation of the results showed that in both the training and validation sets, the risk score was an independent prognostic factor ( Figure 5 ). We then examined the relationship between patients’ risk scores and survival status. In Figures 6A, B , we could find the distribution of patients’ risk scores in ccRCC patients. And, with increasing risk scores, the chance of patient death increased ( Figures 6C, D ). Next, to validate the accuracy of the model, we plotted ROC curves for 1, 2, 3, 4 and 5 years for both datasets. We found the area under the curve (AUC) was found to be almost greater than 0.7 for both datasets from 1 to 5 years, suggesting that the model had good stability and accuracy in predicting patient prognosis ( Figures 6E, F ). Finally, to further test the credibility of the model, we performed a survival analysis in the ccRCC patients ( Figures 7A, B ). At the same time, a further, more specific categorical survival analysis was carried out for all ccRCC patients. The results showed a more rapid decline in survival of ccRCC patients in the high-risk group, irrespective of age, gender grade and stage ( Figures 7C–J ).

As shown in the above analysis, patients in the high-risk group had significantly poorer survival. We therefore wanted to investigate whether there were differences in immune function in order to guide the treatment of the disease in some sense. The results showed more immune cell infiltration in the high-risk group, consisting of T cells, B cells and macrophage cells ( Figure 8A ). Furthermore, almost all immune checkpoint genes were also more highly expressed in the high-risk group ( Figure 8B ), indicating that it is possible that high-risk group ccRCC patients may receive more benefit from immunotherapy.

In order to better predict the prognosis of ccRCC patients, a nomogram was constructed including clinical information and risk score. In Figure 9A , with the use of gender, age, total stage, M stage, grade and risk score values for the patient “TCGA-CZ-4853”, we predicted his mortality rates of 0.0804, 0.207 and 0.325 at 1, 3 and 5 years. Next, we constructed a calibration curve ( Figure 9B ) and found that the nomogram was a good predictor of prognosis at 1, 3 and 5 years for ccRCC patients. In addition, ROC analysis was carried out to better assess the accuracy of the nomogram. The results showed that both the 1 year, 3 year, and 5 years, nomogram was more accurate than clinical information ( Figures 9C–E ).

To assess the knockdown efficiency of the VMP1 in 786-O cells, we examined the expression of the VMP1 in the 786-O cell line by qRT-PCR. Figure 10A showed significant downregulation of VMP1 expression levels in the 786-O cells after siRNA transfection, suggesting that further studies are feasible and meaningful. The EdU assay was used to test whether VMP1 knockdown had an effect on the proliferation of 786-O cells, which showed that the proliferation of cells was suppressed after the vmp1 gene was knocked down ( Figure 10B ).

The effect of decreased VMP1 expression on cell migration and invasion was examined by the transwell method. Results displayed that the decreased expression of VMP1 also impaired cell migration and invasion ( Figures 10C, D ). Scratch healing assays also showed a significantly slower wound healing rate in 786-O cells with a decreased expression of the VMP1 ( Figures 10E, F ).

We analyzed the influence of VMP1 on the apoptosis of 786-O cells. The results indicated that the apoptosis level in the low VMP1 expression group was significantly higher compared to NC group ( Figures 10G, H ).