Differentially expressed autophagy-related genes are potential prognostic and diagnostic biomarkers in clear-cell renal cell carcinoma

Identification of differentially expressed ARGs in ccRCC tissue samples

We analyzed the expression of 232 ARGs in 539 ccRCC and 72 normal kidney tissue samples using the Wilcoxon signed-rank test and identified 38 differentially expressed ARGs (DEARGs). This included 31 upregulated DEARGs and 7 downregulated DEARGs (FDR < 0.05, |log2FC| > 1; Figure 1).

Figure 1. Differential expression of autophagy-related genes in ccRCC tissue samples. The differential expression of 238 autophagy related genes (ARGs) in ccRCC tissue samples (n=539) compared with normal healthy kidney samples (n=72) is shown in the –log (FDR) vs. log (FC) plot. The red dots represent 31 upregulated DEARGs, the green dots represent 7 downregulated DEARGs, and the remaining black dots represent ARGs that are not differentially expressed in ccRCC tissue samples.

Identification of prognostic risk DEARGs in the training group ccRCC patients

Next, we performed univariate Cox regression analysis of the expression of the 38 DEARGs in the training group ccRCC patients to identify prognostic DEARGs. The data showed that expression of 19 DEARGs each significantly correlated with the OS (p < 0.05) and the DFS (p < 0.05) of ccRCC patients (Supplementary Figure 1). We performed Lasso regression analysis to eliminate false-positive DEARGs. Subsequently, our analysis showed that the expression of 12 DEARGs correlated with the OS of ccRCC patients, and the expression of 6 DEARGs correlated with the DFS of ccRCC patients (Supplementary Figure 2). We used these 18 DEARGs (12 for OS and 6 for DFS) to construct models for predicting the prognosis of ccRCC patients.

To determine the optimal model for predicting prognosis, we performed multivariate Cox proportional hazards regression analysis using forward and backward selection algorithms. We identified PRKCQ, BID, BAG1, BIRC5, and ATG16L2 as risk genes in the OS model and EIF4EBP1, BAG1, ATG9B, and BIRC5 as risk genes in the DFS model. The high-risk DEARGs negatively correlated with patient prognosis, whereas the low-risk DEARGs positively correlated with patient prognosis. We identified BID, BIRC5, and ATG16L2 as high risk genes and BAG1 and PRKCQ as low risk genes in the OS model; similarly, we identified EIF4EBP1, ATG9B, and BIRC5 as high risk and BAG1 as a low risk gene in the DFS model (Figure 2A–2B and Table 1).

Figure 2. Characteristics of risk DEARGs in the prognostic risk models. Regression coefficients and hazard ratios of the risk DEARGs for the (A) Overall survival (OS) and (B) Disease-Free Survival (DFS) models are shown.

Testing the prognostic risk models in the training group

We used the following formula to calculate the prognostic risk scores for the training groups using the gene expression values and the regression coefficients of the risk genes: Training group risk score for OS = (-0.1504 × expression value of PRKCQ) + (0.1673 × expression value of BID) + (-0.0944 × expression value of BAG1) + (0.1141 × expression value of BIRC5) + (0.0727 × expression value of ATG16L2); Training group risk score for DFS= (0.0103 × expression value of EIF4EBP1) + (0.1197 × expression value of ATG9B) + (-0.0923 × expression value of BAG1) + (0.2991 × expression value of BIRC5).

Patients in the training group were subdivided into high and low risk groups for OS (n = 133 each) and DFS (n=108 each) based on the median risk scores. Kaplan-Meier survival curve analysis using the log-rank test showed that the high-risk groups exhibited worse prognosis than the low-risk groups (Figure 3A–3B; all p < 0.05). The 3-year and 5-year OS rates for the high-risk patients were 63.3% and 46.0%, respectively, whereas, the 3-year and 5-year OS rates for the low-risk patients were 87.4% and 80.1%, respectively. The 3-year and 5-year DFS rates for the high-risk patients were 58.3.7% and 44.6.8%, respectively. The 3-year and 5-year DFS rates for the low-risk training group patients were 91.0% and 85.6%, respectively.

Figure 3. Analysis of OS and DFS prognostic risk models in training group ccRCC patients. (A) Kaplan-Meier survival curve analysis of OS in the high-risk (red line) and low-risk (green line) ccRCC patients in the training group. (B) Kaplan-Meier survival curve analysis of disease-free survival (DFS) in the high-risk (red line) and low-risk (green line) ccRCC patients. (C) Time-dependent ROC curves show area under curve (AUC) values at 3-year (blue) and 5-year (red) OS in the training group ccRCC patients. (D) Time-dependent ROC curves show AUC values at 3-year (blue) and 5-year (red) DFS in the training group ccRCC patients.

We then measured the predictive performance of the prognostic risk models for the 3-year and 5-year OS and DFS using the time-dependent receiver operating characteristic (ROC) curves. The area under the ROC (AUC) values for the two prognostic models were 0.718 (OS) and 0.779 (DFS) at 3 years, and 0.752 (OS) and 0.783 (DFS) at 5 years (Figure 3C–3D). We then ranked the risk scores of patients for OS and DFS and analyzed their distribution (Figure 4A–4B). The dot plots show the OS and DFS status of individual ccRCC patients in the training groups (Figure 4C–4D). The heat maps show the expression patterns of the risk genes in the high- and low-risk patient groups (Figure 4E–4F). As shown, patients with high-risk scores in the OS group showed upregulation of BID, BIRC5, and ATG16L2 (high risk genes) and downregulation of PRKCQ and BAG1 (protective gene), whereas, patients with low-risk scores showed downregulation of BID, BIRC5, and ATG16L2 and upregulation of PRKCQ and BAG1. Moreover, patients with high-risk groups in the DFS group showed upregulation of EIF4EBP1, ATG9B, and BIRC5 (high risk genes) and downregulation of BAG1 (protective gene), whereas, patients with low-risk groups showed downregulation of EIF4EBP1, ATG9B, and BIRC5 and upregulation of the BAG1.

Figure 4. Prognosis of high-risk and low-risk training group ccRCC patients. (A) Risk score distribution of high-risk (red) and low-risk (green) ccRCC patients in the OS model. (B) Risk score distribution of high-risk (red) and low-risk (green) ccRCC patients in the DFS model. (C) Scatter plot shows the survival status of ccRCC patients in the OS model. Red dots denote patients that are dead and green dots denote patients that are alive. (D) Scatter plot shows survival status of ccRCC patients in the DFS model. Red dots denote patients that are dead and green dots denote patients that are alive. (E) Expression of risk genes in the high-risk (blue) and low-risk (pink) training group ccRCC patients in the OS model. (F) Expression of risk genes in the high-risk (blue) and low-risk (green) training group ccRCC patients in the DFS model. The color code for gene expression in E and F shows green denoting lowest expression and red denoting highest expression

Verification of the prognostic models in the testing group

We validated the accuracy of the two prognostic risk models by analyzing the patients in the testing group. We calculated the risk scores of each patient based on the expression of the risk genes for the OS model, namely, PRKCQ, BID, BAG1, BIRC5, and ATG16L2, and the risk genes for the DFS model, namely, EIF4EBP1, BAG1, ATG9B, and BIRC5. We subdivided the testing group patients into high-risk (n = 130) and low-risk (n = 134) groups in the OS model, and high-risk (n = 113) and low-risk (n = 102) groups in the DFS model. Kaplan-Meier survival curve analysis showed significant differences between the high- and the low-risk groups for OS (p < 0.05) and DFS (p < 0.05), respectively (Figure 5A–5B). Our analysis showed that the 3-year and 5-year OS rates for the high-risk group were 66.8% and 50.2%, respectively, and 89.0% and 81.0%, respectively, for the low-risk groups. The 3-year and 5-year DFS rates for the high-risk group were 65.1% and 58.1%, respectively, and 90.3% and 87.1%, respectively, for the low-risk group. We further performed ROC curve analyses of OS and DFS at 3 years and 5 years for the testing patients based on the two prognostic risk models. The AUC values for the OS model at 3 years and 5 years were 0.695 and 0.709, respectively, whereas, the AUC values for the DFS model at 3 years and 5 years were 0.731 and 0.734, respectively (Figure 5C–5D). The distribution of risk score, survival status, and the expression of risk genes in the testing group patients are shown in Figure 6A–6F.

Figure 5. Validation of the OS and DFS prognostic risk models in the testing group ccRCC patients. (A) Kaplan-Meier survival curve analysis of OS in the high-risk (red line) and low-risk (green line) ccRCC patients in the testing group. (B) Kaplan-Meier survival curve analysis of DFS in the high-risk (red line) and low-risk (green line) ccRCC patients in the testing group. (C) Time-dependent ROC curve analyses shows AUC values for 3-year (blue) and 5-year (red) OS in the testing group ccRCC patients. (D) Time-dependent ROC curve analyses shows AUC values for 3-year (blue) and 5-year (red) DFS in the testing group ccRCC patients.

Figure 6. Prognostic analyses of high-risk and low-risk ccRCC patients in the testing group. (A) Risk score distribution of high risk (red) and low-risk (green) ccRCC patients from the testing group using the OS model. (B) Risk score distribution of high risk (red) and low-risk (green) ccRCC patients from the testing group using the DFS model. (C) Scatter plots show survival status of testing group ccRCC patients using the OS model. (D) Scatter plots show survival status plots of testing group ccRCC patients using the DFS model. (E) Expression of risk genes in the high-risk (blue) and low-risk (pink) testing group ccRCC patients in the OS model. (F) Expression of risk genes in the high-risk (blue) and low-risk (pink) testing group ccRCC patients in the DFS model. The color code for gene expression in E and F shows green denoting lowest expression and red denoting highest expression.

Our data showed that for both OS and DFS, the protective genes were upregulated and the risk genes were downregulated in the low-risk groups, whereas, the protective genes were downregulated and the risk genes were upregulated in the high-risk groups. These results confirmed that both the risk models accurately predicted the prognosis of ccRCC patients

Both prognostic risk models are independently associated with OS and DFS of ccRCC patients from the TCGA group

We further analyzed the relationship between the risk scores of the two models and clinical parameters such as age, gender, histological grade, and the pathological. stage. As shown in Table 2, univariate Cox regression analysis showed that age, histological grade, pathological stage, and the risk score correlated with OS of ccRCC patients (p < 0.05). Similarly, histological grade, pathological stage, and risk score associated with DFS of ccRCC patients (p < 0.05). Next, we performed multivariate Cox regression analysis using age, gender, histological grade, pathological stage, and the risk scores obtained from the prognostic models as explanatory variables. As shown in Table 2, both the prognostic models independently correlated with OS and DFS of all patients from the TCGA group (p < 0.05). Furthermore, multivariate analysis showed that histological grade and pathological stage significantly correlated with OS and DFS (p < 0.05). Moreover, age of the patients significantly correlated with OS (p < 0.05). These results demonstrate that both prognostic models can be independently used to predict OS and DFS in ccRCC patients.

To better predict prognosis at 3- and 5-years post-surgery for ccRCC patients, we constructed two new nomograms using variables associated with OS (age, histological grade, pathological stage, and risk score) and DFS (histological grade, pathological stage, and risk score), respectively (Figure 7A–7B). We then assessed the accuracy of the nomograms using the ROC curve analysis to determine the area under ROC curve (AUC) values. The AUC values for 3-year OS and DFS were 0.812 and 0.871, respectively, whereas the AUC values for 5-year OS and DFS were 0.774 and 0.844, respectively (Figure 7C–7D). These data suggested that both nomograms accurately predicted the 3-year and 5-year OS and DFS rates after surgery in the ccRCC patients.

Figure 7. Construction of nomograms and ROC curve analysis of prognosis for ccRCC patients from the TCGA database. (A–B) The nomograms for (A) OS and (B) DFS are shown. (C–D) ROC curve analysis shows 3-year (blue) and 5-year (red) OS and the corresponding AUC values for the ccRCC patients from the TCGA database. (D) ROC curve analysis shows 3-year (blue) and 5-year (red) DFS and the corresponding AUC values for the ccRCC patients from the TCGA database.

Diagnostic value of risk DEARGs in ccRCC

We then assessed the diagnostic values of the risk genes for OS (PRKCQ, BID, BAG1, BIRC5, and ATG16L2) and DFS (EIF4EBP1, BAG1, ATG9B, and BIRC5) using ROC curve analysis of gene expression data from 535 ccRCC patients from the TCGA database and 28 healthy individuals from the Genotype-Tissue Expression project. The AUC values for PRKCQ,BID, BAG1,BIRC5, ATG16L2, EIF4EBP, and ATG9B genes were 0.649, 0.727, 0.955, 0.868, 0.942, 0.909 and 0.619, respectively (Figure 8). This demonstrated that risk DEARGs are potential diagnostic markers.

Figure 8. ROC curve analysis to determine potential diagnostic value of the risk DEARGs in ccRCC. The ROC curve plots for (A) PRKCQ (AUC = 0.649), (B) BID (AUC = 0.727), (C) BAG1 (AUC = 0.955), (D) BIRC5 (AUC = 0.868), (E) ATG16L2 (AUC = 0.942), (F) EIF4EBP1 (AUC = 0.909), and (G) ATG9B (AUC = 0.619) genes in ccRCC are shown.

Drug sensitivities of ccRCC cell lines

Next, we used the GDSC database to analyze the relationship between the drug sensitivity of several ccRCC cell lines and the expression of BAG1 and BIRC5, which are risk DEARGs for OS and DFS, respectively. Initially, we analyzed the relationship between the expression of these two genes and the IC50 (natural log half maximal inhibitory concentration) values of several targeted drugs in the ccRCC cell lines. We hypothesized that positive correlation between the risk DEARG expression and the IC50 value would imply increased drug resistance in the ccRCC cell lines. Conversely, negative correlation between the risk DEARG expression and IC50 value would suggest increased drug sensitivity in the ccRCC cell lines. Our results indicated that high expression of BAG1 increased the resistance of ccRCC cell lines to drugs such as Bortezomib, Idelalisib, Shikonin, YM201636, Cabozantinib, and NG-25 (p < 0.05), and increased sensitivity of ccRCC cell lines to other drugs such as Erlotinib, AZ628, Lapatinib, A-770041, and HG-5-88-01(p < 0.05; Figure 9A). On the other hand, high expression of BIRC5 expression increased resistance of ccRCC cell lines to drugs such as Salubrinal, PHA-665752, GNF-2, Imatinib, Nilotinib, and Selumetinib (p < 0.05), and increased sensitivity of ccRCC cell lines to other drugs such as CGP-60474, BMS-536924, JW-7-52-1, Panobinostat, Dasatinib, WH-4-023, Rapamycin, Saracatinib, and Bryostatin 1(p < 0.05; Figure 9B).

Figure 9. Correlation between the expression status of risk DEARGs and drug sensitivity of ccRCC cell lines. The plot shows the correlation between the expression status of (A) BAG1 and (B) BIRC5 genes relative to the sensitivity of several ccRCC cell lines to various drugs. The green dots represent drugs that negatively correlate with the expression of the risk genes (p < 0.05) based on their IC50 values; red dots indicate positive correlation of the corresponding drugs with the expression of risk genes (p < 0.05) based on their IC50 values; and black dots represent drugs that do not show any significant correlation based on their IC50 values with the expression of risk genes (p > 0.05).

We also compared the sensitivity of ccRCC cell lines to a BIRC5-targeted drug (Sepantronium bromide) against conventional targeted drugs (Axitinib and Cabozantinib). Our data suggested that the sensitivity of ccRCC cell lines to the BIRC5-targeted drug, Sepantronium bromide, was significantly higher than the sensitivity to conventional targeted drugs (Figure 10A–10C; Table 3; p < 0.05). These results indicate that BIRC5 is a potential therapeutic target for ccRCC.

Figure 10. Drug sensitivity analyses of ccRCC cell lines. The AUC versus IC50 plots show sensitivity of several ccRCC cell lines to treatment with (A) Sepantronium bromide, (B) Axitinib, and (C) Cabozantinib. The cell lines with IC50 values that are greater than the maximum screening concentrations used for the targeted drugs are considered to be resistant to the corresponding drugs. The green dots denote drug-sensitive ccRCC cell lines and red dots denote drug-resistant ccRCC cell lines. IC50 =half maximal inhibitory concentration; AUC: Area under the dose-response curve.

Patient information and databases

We obtained a list of 232 human autophagy-related genes from the Human Autophagy Database (HADb; http://autophagy.lu/clustering/index.html). The gene expression profiles (HTSeq - FPKM) of 611 patients with ccRCC were downloaded from the TCGA database (https://portal.gdc.cancer.gov/), and the patient clinical data was obtained from the cBioPortal database (https://www.cbioportal.org/). The data was analyzed using the R software (https://www.r-project.org/). First, we matched the ID numbers of the patients with their corresponding gene expression profile and clinical data and excluded patients whose ID numbers did not match. We obtained complete gene expression profiles and overall survival (OS) data from 530 patients and complete gene expression profiles and disease-free survival (DFS) data from 431 patients (Table 4). We used the UCSC database (https://xena.ucsc.edu/) to download gene expression profiles (RNA-seq) from the normal kidney tissues of healthy individuals from the Genotype-Tissue Expression project (n=28) and ccRCC tissues from the TCGA database (n = 535). To ensure a unified standard, the RNA-seq profiles were transformed using the formula log₂(x+1) and normalized. The drug sensitivity data for the ccRCC cell lines was downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/).

Identification of DEARGs

We identified the differentially expressed ARGs (DEARGs) using the Wilcoxon signed-rank test. The cut-off values were determined according to the parameters, false discovery rate (FDR) < 0.05 and |log₂ fold change (FC)| > 1.

Construction of OS and DFS risk prognostic models

We randomly divided the 530 patients with complete OS information into two groups, namely, a training group of OS (n = 266) and a testing group of OS (n = 264; Table 5). We also randomly divided the 431 patients with complete DFS information into two groups, namely, a training group of DFS (n = 216) and a testing group of DFS (n = 215; Table 5). We constructed the Cox regression models for OS and DFS using the data from the training group and validated its accuracy using the testing group. Initially, the potential prognostic genes were selected using the univariate Cox regression analysis. Lasso regression analysis was used to eliminate false positives because of over-fitting. Finally, Cox proportional hazards regression was used to construct the OS and DFS prognostic risk models.

Risk score calculation

The risk score for each patient was calculated using the regression coefficients of the individual genes obtained from the multivariate Cox regression model and the expression value of each of the selected genes. The computational formula used for this analysis was

where gene_i represents the identifier of the _ith selected genes and coefficient (gene_i) value is the estimated regression coefficient of the gene_i derived from the Cox proportional hazards regression analysis. The risk score is a measure of prognostic risk for each ccRCC patient. We classified the ccRCC patients into high-risk and low-risk groups using the median risk score of the training group as the cutoff. A high-risk score indicates poor prognosis for the ccRCC patients.

Statistical analysis

All statistical analyses were performed using the R software, and p < 0.05 was regarded as statistically significant. The distribution differences among the variables were assessed using the chi-square test or the Fisher’s exact test. The Kaplan-Meier survival curve analysis and the log-rank test were used to analyze OS and DFS. The Cox regression model was used to analyze the factors that affect the survival of ccRCC patients. Univariate and multivariate analyses were performed using the Cox proportional hazards regression model. Time-dependent ROC analysis was used to evaluate the accuracy of the models that predicted prognosis. ROC curve analysis was also used to estimate the diagnostic value of gene expression. An area under the curve (AUC) value of 0.75 or greater was considered a significant predictive value, and values equal to or greater than 0.6 were regarded as acceptable for predictions.