Prognostic value of immune-related genes in clear cell renal cell carcinoma

Expression of IRGs in ccRCC

The mRNA levels of 2498 IRGs in ccRCC (n = 539) and normal kidney tissues (n = 72) in The Cancer Genome Atlas (TCGA) were examined, and these values were compared through the Wilcoxon signed-rank test. This analysis revealed 681 differentially expressed IRGs (DEIRGs), including 565 genes that were upregulated and 116 genes that were downregulated in ccRCC tissues compared with normal kidney tissues (false-discovery rate [FDR] < 0.05, |log2 fold-change [FC]| > 1) (Figure 1).

Identification of prognostic DEIRGs

To identify possible prognostic DEIRGs (PDEIRGs), we performed a univariate Cox regression analysis of the expression of each DEIRG in the entire TCGA cohort. In total, 263 DEIRGs were found to be significantly associated with the overall survival (OS) of ccRCC patients (p < 0.05) (Figure 2).

Construction of a transcription factor regulatory network

To determine the possible mechanisms behind the dysregulation of PDEIRG expression in ccRCC, we analyzed the correlation between cancer transcription factor (TF) and PDEIRG expression. First, we examined the mRNA levels of TFs in ccRCC (n = 539) and normal kidney tissues (n = 72), and identified 60 TFs (FDR < 0.05, |log2 FC| > 1) that were significantly differentially expressed between the two tissue types (Figure 3A and 3B). Next, we analyzed the correlations between the mRNA levels of the 60 TFs and the PDEIRGs, using a correlation coefficient > 0.4 and a p-value < 0.05 as the cut-off values. Among the 60 TFs, 38 were prominently associated with the aberrant expression of PDEIRGs (p < 0.05). To better explain the regulatory relationships, we constructed a TF-based regulatory network, as displayed in Figure 3C.

Training cohort to identify prognostic genes for inclusion in the risk model

Considering the impact of the PDEIRGs on the OS of patients, we further screened the PDEIRGs to construct a Cox regression hazards model. For this analysis, we used 266 of the 530 patients as a training cohort. First, to avoid overfitting the model, we used Lasso regression to delete PDEIRGs that correlated highly with one another. We thus obtained 12 candidate PDEIRGs (Figure 4A and 4B), and further analyzed them via multivariate Cox proportional hazards regression analysis (with forward selection and backward selection). Ultimately, we obtained seven optimal PDEIRGs (risk genes) for inclusion in the prognostic risk model: PLAU, ISG15, IRF9, ARG2, RNASE2, SEMA3G and UCN. Among these genes, PLAU, ISG15, IRF9, ARG2, RNASE2 and UCN were identified as high-risk genes (predicting a poor prognosis), while SEMA3G was identified as a low-risk gene (serving as a protective factor) in terms of the OS of patients (Figure 5).

Construction of the prognostic risk model in the training cohort

To investigate the significance of the risk genes in predicting the prognosis of ccRCC patients, we used the mRNA levels and estimated regression coefficients of the risk genes to calculate a risk score for each patient. The computational formula was as follows:

Training cohort risk score = (0.0142 × expression of PLAU) + (0.0068 × expression of ISG15) + (0.1467 × expression of IRF9) + (0.0078 × expression of ARG2) + (0.0369 × expression of RNASE2) + (-0.0743 × expression of SEMA3G) + (0.1268 × expression of UCN).

According to the median risk score, the patients in the training cohort were sorted into a high-risk group (n = 133) and a low-risk group (n = 133). To determine the prognostic difference between the high-risk and low-risk groups, we created a Kaplan-Meier curve based on the log-rank test. The prognosis was worse in the high-risk group than in the low-risk group (p < 0.05) (Figure 6A). The OS rates at three years and five years for the high-risk group in the training cohort were 63.1% and 41.4%, respectively, while the corresponding rates for the low-risk group were 90.8% and 87.1%, respectively. We then used time-dependent receiver operating characteristic (ROC) curves to examine the predictive accuracy of the model for OS at three years and five years. The area under the ROC (AUC) values for the prognostic model were 0.760 at three years and 0.789 at five years (Figure 6B). We then ranked the risk scores of the patients in the training cohort and analyzed their distribution (Figure 6C). The survival status of each patient in the training cohort is marked on the dot plot in Figure 6D.

A heat map was generated to depict the expression patterns of the risk genes in the two prognostic groups (Figure 6E). In patients with high risk scores in the training cohort, the six high-risk genes (PLAU, ISG15, IRF9, ARG2, RNASE2 and UCN) were upregulated, while the protective gene (SEMA3G) was downregulated. In patients with low risk scores, these risk genes displayed the opposite expression pattern.

Verification of the performance of the prognostic model

To verify the accuracy of the prognostic risk model, we used it to analyze the testing cohort (the remaining 264 patients from the 530 total) and the entire TCGA cohort. First, we used the seven risk genes (PLAU, ISG15, IRF9, ARG2, RNASE2, SEMA3G and UCN) to calculate the risk score of each patient in the testing cohort and the entire TCGA cohort. The patients in each cohort were then classified into two groups based on the median risk score of the training cohort. In the testing cohort, 148 patients were categorized as high-risk and 116 were categorized as low-risk. In the entire TCGA cohort, 290 patients were classified as high-risk and 240 were classified as low-risk.

Next, Kaplan-Meier survival analysis was used to determine the prognostic differences between the high-risk and low-risk groups. The Kaplan-Meier survival curves differed significantly between the two risk groups in both the testing cohort and the entire TCGA cohort (p < 0.05) (Figure 7A and 7B); throughout the follow-up time, the survival rate was higher in low-risk patients than in high-risk patients. In the testing cohort, the survival rates at three and five years in the high-risk group were 66.2% and 51.7%, respectively, while the corresponding rates in the low-risk group were 89.0% and 81.0%, respectively. In the entire TCGA cohort, the survival rates at three and five years in the high-risk group were 66.3% and 50.8%, respectively, while the corresponding rates in the low-risk group were 89.2% and 80.6%, respectively. ROC analyses were performed for the testing cohort and the entire TCGA cohort at three and five years. In the testing cohort, the AUCs at three and five years were 0.715 and 0.692, respectively. In the entire TCGA cohort, the AUCs at three and five years were 0.740 and 0.745, respectively (Figure 7C and 7D).

The risk score distribution, survival status and risk gene expression in the testing cohort and the entire TCGA cohort are displayed in Figure 8A–8F. Similar to the results in the training cohort, protective gene levels were higher and risk gene levels were lower in the low-risk group than in the high-risk group. These results indicated that our prognostic risk model is capable of precisely predicting the prognosis of ccRCC patients.

Independent prognostic value of the risk model in the entire TCGA cohort

Next, we performed univariate and multivariate Cox regression analyses to assess whether the risk score generated by our model was independent from other clinical parameters (age, gender, histological grade and pathological stage) as a prognostic factor for ccRCC. The univariate analysis indicated that the variables of age, histological grade, pathological stage and risk score were associated with the prognosis of ccRCC patients. The multivariate analysis revealed that the risk score was independently associated with OS in the entire TCGA cohort (p < 0.05) (Table 1). These results demonstrated that the prognostic risk model can be used independently to predict the prognosis of ccRCC patients. However, the three clinical variables (age, histological grade and pathological stage) were also found to be significant prognostic factors in the multivariate analysis (p < 0.05).

We then assessed whether the risk score from our model was more accurate than the other clinical parameters (age, histological grade and pathological stage) in predicting OS at three and five years. Indeed, the risk score was more accurate than the other clinical parameters: the AUCs at three years for age, histological grade and pathological stage were 0.568, 0.676 and 0.689, respectively (Figure 9A), and the corresponding values at five years were 0.587, 0.642 and 0.643, respectively (Figure 9B).

To better predict the prognosis of ccRCC patients at three and five years post-surgery, we constructed a new nomogram from the variables associated with OS (age, histological grade, pathological stage and risk score) (Figure 10A). ROC curve analysis was used to evaluate the accuracy of the nomogram. The nomogram could accurately predict OS at three and five years post-surgery, with AUCs of 0.814 and 0.775, respectively (Figure 10B).

Clinical utility of the prognostic risk model

To examine the ability of our model to predict the progression of ccRCC, we analyzed the relationships between the risk factors from our model (the risk score and risk genes) and the clinical variables (age, gender, histological grade and pathological stage) in the entire TCGA cohort. As the values of certain factors increased (PLAU, RNASE2 and UCN levels and the risk score), the histological grade of ccRCC patients increased (all p < 0.05) (Figure 11A–11D), and as the values of other factors increased (PLAU, ISG15, IRF9, RNASE2 and UCN levels and the risk score), the pathological stage of ccRCC patients increased (all p < 0.05) (Figure 11E–11J). The expression of UCN was higher in patients > 60 years old than in those ≤ 60 years old (p < 0.05) (Figure 11K). In contrast, as SEMA3G expression increased, the values of two clinical variables (histological grade and pathological stage) decreased (both p < 0.05) (Figure 11L and 11M) (Table 2). These results demonstrated that the dysregulation of immune-related risk gene expression is associated with the development of ccRCC.

To determine whether our model could reflect the status of the tumor immune microenvironment in patients, we analyzed the correlation between the risk score and immune cell infiltration in the entire TCGA cohort. As the risk score increased, the content of immune cells (CD8+ T cells, neutrophils, macrophages and dendritic cells) in ccRCC tissues also increased (p < 0.05) (Figure 12A–12F).

Previous studies have indicated that the tumor mutation burden is significantly associated with the clinical effectiveness of immunotherapy [12]. To assess whether our model could stratify patients with different sensitivities to immunotherapy, we compared the mutation counts of the high-risk and low-risk groups in the entire TCGA cohort. The mutation count was greater in the high-risk group than in the low-risk group (p < 0.05) (Figure 13).

Databases

The 2498 IRGs were obtained from the ImmPort database (https://www.immport.org/home). The 539 ccRCC patients’ data, including transcriptomic data, mutation data and clinical information, were downloaded from the TCGA portal (https://portal.gdc.cancer.gov/) and the cBio Cancer Genomics portal (https://www.cbioportal.org/) [28, 29]. All data were processed with R software (https://www.r-project.org/). We matched patients’ transcriptomic data and clinical information according to their ID numbers, and we removed patients if their ID numbers did not match. We thus obtained data from 530 patients with complete gene expression profiles and OS information (Table 3). Immune infiltrate data from the ccRCC patients was obtained from the Cistrome project (http://www.cistrome.org/) [30], which contains the abundances of six types of tumor-infiltrating immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells). Data on TFs associated with cancer were also collected from the Cistrome project.

Identification of DEIRGs

The Wilcoxon signed-rank test was used to screen DEIRGs based on the following cut-off values: FDR < 0.05 and |log2 FC| > 1.

Experimental model construction

The 530 patients were randomly divided into two cohorts: a training cohort (n = 266) and a testing cohort (n = 264) (Table 4). The training cohort was used to construct the Cox regression hazards model, and the testing cohort was used to test the performance of the model. Initially, univariate Cox analysis was used to identify possible PDEIRGs. Next, Lasso regression was used to select potential risk genes and eliminate genes that would overfit the model. Finally, we used Cox proportional hazards regression to construct a prognostic risk model.

Risk score calculation

To calculate the risk score for each patient, we used the regression coefficients from the multivariate Cox regression model to weight the expression values of the selected genes. The following computational formula was used for this analysis:

$Riskscore(patient)={\displaystyle \sum }_{i=1}^{n}coefficient(genei)expressionvalueof(genei)$

Here, ‘gene_i’ is the i^th selected gene, and ‘coefficient (gene_i)’ is the estimated regression coefficient of gene_i from the Cox proportional hazards regression analysis. The risk model was used to measure the prognostic risk of each patient with ccRCC. The median risk score of the training cohort was used as the cut-off value to divide all the ccRCC patients into two groups: the high-risk group and the low-risk group. A high risk score indicates a poor prognosis for ccRCC patients.

Statistical analyses

R software was used to perform all statistical analyses, and p < 0.05 was considered statistically significant. The rank correlation among the different variables was assessed with the Pearson correlation coefficient test. Differences between variables were assessed with independent t-tests. Kaplan-Meier curves and log-rank tests were used to analyze the survival data, and univariate Cox regression analysis was used to identify factors affecting the survival of patients diagnosed with ccRCC. Multivariate Cox regression analysis was used to identify independent prognostic factors. Time-dependent ROC analysis was used to evaluate the accuracy of the prognostic prediction model. An AUC > 0.60 was regarded as acceptable for predictions, and an AUC > 0.75 was deemed to have excellent predictive value [31, 32].