Identification of a 6-Gene Signature Associated with Resistance to Tyrosine Kinase Inhibitors: Prognosis for Clear Cell Renal Cell Carcinoma

Background

Renal cell carcinoma (RCC) is the most common urinary system tumor, accounting for 3.7% of all new malignant tumors [1]. In the past few decades, the development of imaging technology has facilitated diagnosis at earlier stages of the disease, thus resulting in significantly prolonged survival times [2,3]. Unfortunately, approximately one-third of patients diagnosed with clear cell RCC (ccRCC) have local or distant metastases at the time of diagnosis [2,4]; 30% to 40% of patients with early local metastases experience relapse or further metastatic disease, even after radical nephrectomy [5,6]. As such, the prognoses of patients with ccRCC associated with metastatic disease remains poor, with a 5-year survival rate of only 12% [2].

Clear cell is the most common subtype of RCC and accounts for approximately three-quarters of all cases [7]. These tumors are highly angiogenic and are frequently associated with von Hippel-Lindau (VHL) gene mutations. Inactivation of the VHL gene increases hypoxia-inducible factor activity, eventually leading to overexpression of vascular endothelial growth factor (VEGF) and platelet-derived growth factor [8,9]. Given the VHL protein’s role in the pathogenesis of tumor cells, metastatic disease treatment now focuses on targeted therapy based on the VEGF-tyrosine kinase inhibitor (TKI); the receptor TKIs sunitinib and sorafenib are currently first-line treatments for this condition [10–12].

While the administration of TKIs has prolonged the median survival time for patients with advanced ccRCC [13–15], almost all patients develop drug resistance after 6 to 15 months of treatment [16]. The molecular mechanisms underlying drug resistance observed in patients with advanced ccRCC remains unclear, and no reliable biomarkers of this condition have been developed [17]. Therefore, there is an urgent need to explore new prognostic models that might predict survival in patients with advanced ccRCC.

Our study identified several differentially expressed genes (DEGs) associated with sunitinib/sorafenib resistance from published datasets that were identified in the Gene Expression Omnibus (GEO) database. We then evaluated the prognosis with regard to each of the DEGs, using findings maintained in The Cancer Genome Atlas (TCGA) database, and constructed a 6-DEGs risk signature. This risk signature can effectively assess the prognoses of patients diagnosed with ccRCC and is particularly accurate for those with metastatic disease.

Material and Methods

IDENTIFICATION OF DEGS ASSOCIATED WITH DRUG RESISTANCE:

mRNA expression profiles in the GSE64052 file were downloaded from the GEO data repository (https://www.ncbi.nlm.nih.gov/geo/). The data file included 14 transplanted tumor samples that were resistant to TKIs (sunitinib or sorafenib) and 14 untreated control samples. The probe of the raw data file (Series Matrix.txt) was annotated with official gene symbols by the GPL570 platform file. The R package “limma” [18] was used for background correction and difference analysis (the cutoff criterion at a false discovery rate [FDR] <0.05 and the absolute value of the log of the fold change [|logFC|] >1). We identified 91 DEGs related to TKI resistance. Metascape was used for further functional enrichment analysis of DEGs [19].

PREPROCESSING OF THE TCGA DATA:

RNA-Seq data (HT-Seq-Count) and associated clinical data were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). Official gene symbols were processed by the annotation file (Homo_sapiens.GRCh38.99.chr.gtf) within the Genecode database (https://www.gencodegenes.org/) [20]. The R package “edge R” algorithm [21] was used to perform background correction and normalization, and the R package “sva” was used to eliminate batch effects and other off-target variations associated with TCGA and GSE64052 data. The final analysis included 79 DEGs, identified in both TCGA and GEO expression matrices, and 463 patients with complete clinical follow-up information (≥90 days), along with complete clinicopathological information.

ESTABLISHMENT OF A PROGNOSTIC MODEL ASSOCIATED WITH TKI RESISTANCE:

The ccRCC patients were randomly divided into a training cohort (n=232) and a validation cohort (n=231) using the R package “caret”. Table 1 shows the clinical information associated with the training cohort, the validation cohort, and the entire cohort. Univariate Cox proportional hazards regression was analyzed with the R package “survival” to identify relationships between the TKI-resistance-associated DEGs with overall survival (OS); those with P<0.05 were selected as candidate variables. Subsequently, the least absolute shrinkage and selection operator (Lasso) Cox regression method analysis, which is an algorithm based on L1-penalized linear regression to prevent over-fitting of the model, was used to screen candidate variables; the coefficients of each DEG were calculated using the R package “glmnet” [22]. The best prognostic markers among the DEGs were established using these methods. The findings were confirmed in the validation cohort and the entire cohort. The risk score formula, based on 6 DEGs, was constructed as follows:

In this formula, n represents the number of DEGs, βi represents the coefficient of each DEG, and Exp_i represents the level of expression of each DEG. Risk scores were calculated for each patient; the cases were then divided into high-risk and low-risk groups according to the median risk score.

EVALUATION OF THE 6-DEGS PROGNOSIS MODEL:

For survival analysis, a Kaplan-Meier survival curve was plotted, and a Wilcoxon test was adopted to assess significant differences between high-risk and low-risk patient groups. Furthermore, the receiver operating characteristic (ROC) curve, the area under ROC (AUC), and survival-status scatter plot were drawn to evaluate the prognosis model’s accuracy in the training cohort, the validation cohort, and the entire cohort.

IDENTIFICATION OF INDEPENDENT PROGNOSTIC FACTORS ASSOCIATED WITH SURVIVAL:

Stratified analysis and Mann-Whitney tests were conducted to identify the discriminatory ability of risk scores concerning various clinical characteristics (i.e., age, living status, sex, grade, clinical stage, and TNM stage). Additionally, univariate and multivariate Cox regression analyses were performed using R package “survival” to assess the risk scores’ prognostic value and associated clinicopathological parameters. Multifactor ROC curves verified the parameters’ accuracy and specificity.

EVALUATION AND VERIFICATION OF NOMOGRAM:

The R package “rms” was used to construct a nomogram that included risk scores and clinicopathological features to predict the progress and prognosis of ccRCC patients at 1, 3, and 5 years. The Harrell consistency index (C-index) was used to test the accuracy of the nomogram; a calibration curve was constructed to test the consistency based on 1-, 3-, and 5-year survival predictions.

STATISTICS AND PLOTTING:

Statistics and plotting were performed using GraphPad Prism 8.2.1 (GraphPad Software Inc., San Diego, CA) and R environment (version 3.6.2, www.r-project.org/).

Results

CONSTRUCTING A 6-DEGS RISK SCORE ASSOCIATED WITH TKI RESISTANCE:

To identify key genes related to TKI resistance, we downloaded raw data associated with 28 transplanted tumor samples from the GEO database and annotated 21 655 genes using the GPL570 platform document. We identified 91 DEGs in a comparison between gene expression levels in 14 transplanted tumor samples, reported as resistant to sorafenib/sunitinib, and those in 14 untreated samples; 29 DEGs were upregulated, and 62 DEGs were downregulated in the TKI-resistant samples. We identified these genes as associated with drug resistance (Supplementary Table 1). Figure 1A shows a volcano plot of these DEGs.

We next used Metascape to perform GO enrichment analysis on these DEGs. For molecular function, these genes mainly enriched in sodium-phosphate symporter activity, extracellular matrix binding, calcium ion binding, and virus receptor activity. Meanwhile, for biological processes, the DEGs were mainly enriched in drug transport, blood vessel development, hormone metabolic process, regulation of anion transport, and other biological processes (Supplementary Figure 1, Supplementary Table 2). To assess the DEGs’ prognostic value in the cases included in the TCGA database, ccRCC patients identified in this cohort were randomly and evenly divided into a training cohort (n=232) and a validation cohort (n=231). For the training cohort, we used Cox proportional hazards regression model to identify 32 genes that were significantly associated with OS (Supplementary Table 3, P<0.05). Furthermore, to establish a general indicator for patient prognosis, we performed Lasso regression analysis on these 32 genes, determined the best penalty parameters through 10 rounds of cross-validation, and finally obtained 6 DEGs significantly correlated with ccRCC (Figure 1B, 1C). Remarkably, the Spearman correlation analysis confirmed the low correlation among the 6 DEGs (Figure 1D), which underscored the effectiveness of this model.

Among these 6 genes, matrix metalloproteinase (MMP)24, angiotensin-converting enzyme (ACE)2, and solute carrier family 44 member 4 (SLC44A4) were associated with a positive prognosis (hazard ratio [HR] <1), while complement component C1R, chromosome 1 open reading frame (C1ORF)-194, and a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS)-15 were associated with a negative prognosis (HR >1) (Figure 1E). A weighted prognostic risk score formula was established, based on the regression coefficient and expression levels of each of the 6 genes; risk score= (−0.1841×MMP24) +(−0.1096×ACE2) +(−0.1030×SLC44A4) +(0.1490×C1R) +(0.2328×C1ORF194) +(0.2841×ADAMTS15). The ccRCC patients were then divided into a high-risk group (n=131) and a low-risk group (n=130) based on the median risk score of 0.945.

EVALUATION AND VERIFICATION OF PROGNOSTIC MODELS:

The Kaplan-Meier log-rank test was used to evaluate the risk score’s accuracy for predicting survival. The results revealed that patients identified as high-risk via this scoring system had lower OS than patients identified as low-risk (Figure 2A). Additionally, the time-dependent ROC curve showed that the AUCs associated with the 6-DEGs signatures at 1, 3, and 5 years were 0.785, 0.756, and 0.748, respectively (Figure 2B). We also constructed a survival-status scatter plot to assess the distribution of risk scores among patients who were alive and those who had died at the time of the study. The results demonstrated that, compared with patients who remained alive, the distribution of scores among those patients who were deceased was substantially more concentrated; these latter patients were mainly associated with high-risk scores (Figure 2C). Similar results were obtained when this analysis was applied to the validation cohort and the overall cohort (Figure 2D–2I). Table 2 shows the specific results from the survival analysis. In summary, the 6-DEGs-based risk score accurately predicted survival outcomes in patients diagnosed with ccRCC.

Since TKI resistance mainly occurs in patients with advanced metastatic disease, we analyzed the outcomes associated with the 6-DEGs risk score in patients with metastatic disease (M1) identified in the entire cohort (n=67). The results revealed that high-risk M1 patients had significantly lower OS and shorter median survival times than low-risk patients (Figure 2J). The 1-, 3-, and 5-year AUCs of M1 patients were 0.705, 0.749, and 0.796, respectively (Figure 2K). Notably, the survival-status scatter plots showed that M1 patients who were deceased were more concentrated within the high-risk score group (Figure 2L). Taken together, the risk model has superior sensitivity and specificity for ccRCC patients with metastatic disease.

STRATIFIED ANALYSIS OF RISK SCORES:

As a further confirmation of the clinical value of the 6-DEGs risk score model, we performed a stratified analysis with patients divided into different subgroups based on demographics and clinical characteristics, including age (<60 or ≥60), sex (female vs. male), survival (living vs. deceased), pathological grade (I+II vs. III+IV), clinical stage (I+II vs. III+IV), T stage (T1+T2 vs. T3+T4), N stage (N0 vs. N1), and M stage (M0 vs. M1). A Mann-Whitney test revealed that deceased patients, those with high pathological grades (III, IV), high clinical stages (III, IV), and high TNM diagnoses (T3, T4, M1, N1) were assigned higher risk scores. Interestingly, there were no significant differences concerning age or sex (Figure 3).

THE 6-DEGS RISK SCORE AS AN INDEPENDENT PROGNOSTIC FACTOR:

We performed univariate Cox regression and multivariate Cox analysis to evaluate the relationships between the 6-DEGs risk score and each demographic and clinicopathological parameter with OS. The results revealed that age, grade, clinical stage, and risk scores were associated with the prognoses of ccRCC patients (P<0.001; Figure 4A). Therefore, we included these 4 factors for further multivariate analysis and found that the covariates age, stage, and risk score were still significant in the multivariate Cox analysis (P<0.001; Figure 4B). Also, using a multifactor ROC curve for comparisons with age, grade, and stage, we found that the accuracy of the 6-DEGs risk for predicting OS was stable throughout the 1-, 3-, and 5-year periods (Figure 4C–4E).

Moreover, we used univariate Cox analysis to compare various factors (age, grade, T stage, and risk score) that can assess the prognosis of metastatic patients (M1). We found that only risk score was significant in a univariate analysis (P<0.001; Figure 5A). The 1-, 3-, and 5-year ROC curves also showed that the risk score had better accuracy than other variables (Figure 5B–5D). Hence, we conclude that the 6-DEGs signature can function as an independent predictor of prognosis in patients diagnosed with ccRCC, especially with metastatic ccRCC.

NOMOGRAM BASED ON 6-DEGS SIGNATURE FOR PREDICTING THE SURVIVAL TIME OF CCRCC PATIENTS:

Given the significance of the risk score and clinical-pathological parameters for predicting survival, we constructed a nomogram that combined the risk score and specific clinical parameters (age, grade, clinical stage) to predict the survival of patients with maximum efficacy (Figure 6A). The results revealed that the prognostic nomogram could accurately estimate the OS at 1, 3, and 5 years in patients diagnosed with ccRCC (Figure 6B–6D). The C-index of the nomogram was 0.784 (95% confidence interval, 0.745–0.823), a result that indicates that the nomogram can distinguish between patients with good versus poor prognoses.

Discussion

The mechanisms underlying TKI resistance may be related to the development of proangiogenic pathways, the tumor microenvironment, epithelial-mesenchymal transition (EMT), single-nucleotide polymorphisms, and/or other induced genetic alterations [17]. These extensive genetic changes are an important foundation affecting drug resistance. Hence, the identification of key biomarkers related to TKI resistance may help researchers to further clarify the mechanisms underlying ccRCC. Although widely used, the current TNM system and related clinicopathological features are ineffective for determining the prognoses of patients with advanced ccRCC [23].

In the current study, we first identified 32 genes associated with TKI resistance via analysis of differential expression in association with patient prognosis. We then introduced Lasso Cox regression analysis and selected 6 genes with good prognostic value and low correlation to establish a risk signature. Encouragingly, Kaplan-Meier survival, ROC, and univariate and multivariate Cox regression analyses verified the accuracy and specificity of the 6-DEGs signature and its applicability in the overall ccRCC cohort. In particular, survival analysis and univariate Cox analysis showed that the 6-DEGs risk score had good accuracy in predicting patients with metastatic disease. Finally, we developed a nomogram that integrated clinical prognostic variables and risk score to provide clinicians with a tool that can more effectively predict the survival of patients diagnosed with ccRCC, particularly those with advanced metastatic disease.

Among the 6 genes selected for the risk score, ACE2 is a homolog of ACE and plays a key role in the renin-angiotensin system [24]. Qian et al. [25] demonstrated that ACE2 overexpression resulted in upregulated expression of E-cadherin and downregulated expression of vimentin. These results indicated that ACE2 could inhibit EMT and reduce the metastatic potential of lung cancer cells. Likewise, Zhang et al. [26] reported that ACE2 promoted the downregulated expression of VEGF-A in breast cancer cells and thereby reduced angiogenesis. Furthermore, ACE2 has high catalytic efficiency and can hydrolyze angiotensin II into Ang 1–7 [27], a protein that inhibits angiogenesis, invasion, and metastasis of tumor cells [28–31]. MMP24 is a member of the membrane-type MMP family of zinc-dependent endopeptidases that act primarily to promote degradation of extracellular matrix [32]. Previous studies revealed the complex biological characteristics of MMPs. While most MMPs promote tumor cell migration, EMT, adhesion, and angiogenesis, others mediate tumor-suppressive effects [32–35]. Sugimoto et al. [36] found that MMP24 may be upregulated in the extracellular matrix of breast cancer cells, thereby promoting tumor invasiveness. SLC44A4 is a member of the solute carrier protein family; the group including SLC44A1–5 are also known as choline transporter-like proteins (CTLs)1–5 [37]. SLC44A4 promotes the synthesis and transport of acetylcholine [38] and the absorption of thiamine pyrophosphate, a phosphorylated form of vitamin B1 [39]. SLC44A4 expression varies significantly in different tumors and different locations. Song et al. [40] demonstrated that inhibiting SLC44A4 expression results in reduced secretion of acetylcholine, which inhibits the growth of lung cancer cells. C1R is a member of the S1 protein family of peptidases; the gene encodes a proteolytic subunit in the C1 complex that contributes to the classical activation pathway of the complement system [41]. In the tumor microenvironment, complement activation could enhance tumor growth and accelerate metastasis [42]. Riihilä et al. [43] reported elevated levels of C1R expression in squamous cell carcinoma of the skin. Inhibition of C1r or C1s in squamous cell carcinoma cells inhibited the activation of both extracellular signal-related kinase 1/2 and Akt; these actions led to reduced tumor growth and angiogenesis in vivo. C1ORF194 is a protein-coding gene. Previous studies demonstrated that a mutant form of C1ORF194 protein disrupted signaling pathways involving Ca2+ homeostasis, ultimately resulting in Charcot-Marie-Tooth disease [44]. However, the role of C1ORF194 in tumors has not been reported. Finally, ADAMTS15 is a multi-domain matrix-associated zinc metalloendopeptidase; Kelwick et al. [45] reported that ADAMTS15 inhibited breast cancer cell metastasis.

Despite the excellent performance of the 6-DGEs signature as a prognostic indicator in cases of advanced ccRCC, it has several limitations and disadvantages. First, because we had no access to treatment data related to drug resistance, the genes were identified primarily based on existing experimental animal data. As such, additional clinical trials are needed to confirm the validity of the 6-DEGs signature in patients with clinical TKI resistance. Second, although the patient cohorts were grouped randomly to reduce bias, the data in the TCGA database were obtained primarily from white patients. It is necessary to verify our findings in patient cohorts featuring other races and ethnicities. Third, while 5 of the 6 genes contributed to the process of tumorigenesis in previous reports, the reported functions of MMP24, SLC44A4, and ADAMTS15 are inconsistent with our findings. Heterogeneity among tumors is commonplace, and the genes associated with different tumors can vary substantially [46,47]; as such, these results might be expected. However, we have not yet examined the functional mechanisms underlying the 6-DEGs signature genes, notably concerning their roles in promoting or inhibiting renal cancer. Additional experiments are needed to address questions associated with mechanisms underlying the roles of the 6 DEGs in mediating TKI resistance in association with ccRCC.

Conclusions

In summary, we used information available in the GEO and TCGA databases to establish a 6-DEGs signature associated with TKI resistance as a potential prognostic indicator for patients with ccRCC. Univariate analysis, multivariate analysis, and nomogram calibration curves supported the strong predictive value of the 6-DEGs-based risk score. As such, the mechanisms underlying differential expression of these 6 genes and their relationship to the pathogenesis of ccRCC should be explored in future studies.